multifield_plasticity.cpp

Go to the documentation of this file.

/**
 * \file multifield_plasticity.cpp
 *
 * Multifield plasticity with contact
 *
 */
 
/* 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/>. */
 
#include <MoFEM.hpp>
 
#include <chrono> //for debugging
#include <boost/tokenizer.hpp>
#include <IntegrationRules.hpp>
#include <BlockMatData.hpp>
#include <MultifieldPlasticity.hpp>
 
struct MeasureTime {
  int time;
  chrono::high_resolution_clock::time_point start;
  chrono::high_resolution_clock::time_point stop;
  MeasureTime() { start = chrono::high_resolution_clock::now(); }
  ~MeasureTime() {
    stop = chrono::high_resolution_clock::now();
    auto duration = chrono::duration_cast<chrono::microseconds>(stop - start);
    cout << "Time taken by function: " << duration.count() << " us." << endl;
  }
};
 
using namespace MoFEM;
using namespace FTensor;
 
using OpLogStrainMatrixLhs = FormsIntegrators<DomainEleOp>::Assembly<
    USER_ASSEMBLE>::BiLinearForm<GAUSS>::OpGradTensorGrad<1, 3, 3, 1>;
using OpStiffnessMatrixLhs = FormsIntegrators<DomainEleOp>::Assembly<
    USER_ASSEMBLE>::BiLinearForm<GAUSS>::OpGradSymTensorGrad<1, 3, 3, 0>;
using OpBodyForce = FormsIntegrators<DomainEleOp>::Assembly<
    USER_ASSEMBLE>::LinearForm<GAUSS>::OpSource<1, 3>;
 
using OpMixDivULhs = FormsIntegrators<DomainEleOp>::Assembly<
    USER_ASSEMBLE>::BiLinearForm<GAUSS>::OpMixDivTimesVec<3>;
using OpLambdaGraULhs = FormsIntegrators<DomainEleOp>::Assembly<
    USER_ASSEMBLE>::BiLinearForm<GAUSS>::OpMixTensorTimesGrad<3>;
using OpMixDivURhs = FormsIntegrators<DomainEleOp>::Assembly<
    USER_ASSEMBLE>::LinearForm<GAUSS>::OpMixDivTimesU<3, 3, 3>;
using OpMiXLambdaGradURhs = FormsIntegrators<DomainEleOp>::Assembly<
    USER_ASSEMBLE>::LinearForm<GAUSS>::OpMixTensorTimesGradU<3>;
using OpMixUTimesDivLambdaRhs = FormsIntegrators<DomainEleOp>::Assembly<
    USER_ASSEMBLE>::LinearForm<GAUSS>::OpMixVecTimesDivLambda<3>;
using OpMixUTimesLambdaRhs = FormsIntegrators<DomainEleOp>::Assembly<
    USER_ASSEMBLE>::LinearForm<GAUSS>::OpGradTimesTensor<1, 3, 3>;
 
using OpSpringLhs = FormsIntegrators<BoundaryEleOp>::Assembly<
    USER_ASSEMBLE>::BiLinearForm<GAUSS>::OpMass<1, 3>;
using OpSpringLhsNoFS = FormsIntegrators<BoundaryEleOp>::Assembly<
    USER_ASSEMBLE>::BiLinearForm<GAUSS>::OpMass<1, 3>;
using OpSpringRhs = FormsIntegrators<BoundaryEleOp>::Assembly<
    USER_ASSEMBLE>::LinearForm<GAUSS>::OpBaseTimesVector<1, 3, 1>;
 
using OpInertiaForce = FormsIntegrators<DomainEleOp>::Assembly<
    USER_ASSEMBLE>::LinearForm<GAUSS>::OpBaseTimesVector<1, 3, 1>;
 
using OpCentrifugalForce2 = FormsIntegrators<DomainEleOp>::Assembly<
    USER_ASSEMBLE>::LinearForm<GAUSS>::OpSource<1, 3>;
 
// using OpScaleL2 = MoFEM::OpScaleBaseBySpaceInverseOfMeasure<DomainEleOp>;
 
#include <BasicFiniteElements.hpp>
#include <AnalyticalSurfaces.hpp>
#include <RigidBodies.hpp>
#include <BasicFeTools.hpp>
#include <MatrixFunction.hpp>
 
// global variables
BlockParamData *cache;
std::map<int, BlockParamData> mat_blocks;
boost::weak_ptr<MoFEM::BlockMatContainer> MoFEM::block_mat_container_ptr;
 
#include <ElasticOperators.hpp>
#include <PlasticOperators.hpp>
#include <ContactOperators.hpp>
#include <RotatingFrameOperators.hpp>
#include <DualBase.hpp>
 
using namespace OpElasticTools;
using namespace OpPlasticTools;
using namespace OpContactTools;
using namespace OpRotatingFrameTools;
 
// for testing
constexpr bool TEST_H1_SPACE = false;
EntityType zero_type = TEST_H1_SPACE ? MBVERTEX : MBTET;
FieldSpace space_test = (TEST_H1_SPACE ? H1 : L2);
 
struct ContactPlasticity {
 
  ContactPlasticity(MoFEM::Interface &m_field) : mField(m_field) {}
 
  MoFEMErrorCode runProblem();
  MoFEMErrorCode loadMeshBlockSets();
  static MoFEMErrorCode getAnalysisParameters();
  static MoFEMErrorCode myMeshPartition(MoFEM::Interface &m_field);
 
private:
  MoFEM::Interface &mField;
  MoFEMErrorCode setUP();
  MoFEMErrorCode createCommonData();
  MoFEMErrorCode OPs();
  MoFEMErrorCode bC();
  MoFEMErrorCode tsSolve();
  MoFEMErrorCode postProcess();
  MoFEMErrorCode checkResults();
 
  MatrixDouble invJac, jAc;
  boost::shared_ptr<OpContactTools::CommonData> commonDataPtr;
  boost::shared_ptr<PostProcVolumeOnRefinedMesh> postProcFe;
  std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uXScatter;
  std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uYScatter;
  std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uZScatter;
 
  // mofem boundary conditions
  boost::ptr_map<std::string, NeumannForcesSurface> neumann_forces;
  boost::ptr_map<std::string, NodalForce> nodal_forces;
  boost::ptr_map<std::string, EdgeForce> edge_forces;
  boost::shared_ptr<DirichletDisplacementRemoveDofsBc> dirichlet_bc_ptr;
 
  boost::shared_ptr<std::vector<unsigned char>> boundaryMarker;
  boost::shared_ptr<Range> allBoundaryEnts;
  // for forces scaling
  boost::ptr_vector<DataFromMove> bcData;
  boost::shared_ptr<DomainSideEle> volSideFe;
 
  // boost::shared_ptr<DomainEle> contactVolFeLhs;
  // boost::shared_ptr<DomainEle> contactVolFeRhs;
  boost::shared_ptr<DomainEle> feLambdaRhs;
  boost::shared_ptr<DomainEle> feLambdaLhs;
 
  Range getEntsOnMeshSkin();
 
  // SmartPetscObj<DM> dmElastic;
  SmartPetscObj<DM> dmEP;
  SmartPetscObj<DM> dmTAU;
  SmartPetscObj<DM> dmContact;
 
public:
  // global params
  struct ProblemData {
    int order;
    PetscBool no_output;
    PetscBool debug_info;
    int output_freq;
    int atom_test_nb;
    int reaction_id;
    PetscBool scale_params;
 
    PetscBool is_large_strain;
    PetscBool is_fieldsplit_bc_only;
 
    PetscBool is_neohooke;
    PetscBool with_plasticity;
    PetscBool calculate_reactions;
    PetscBool is_ho_geometry;
    PetscBool save_restarts;
    PetscBool is_restart;
 
    PetscBool is_lambda_split;
 
    PetscBool with_contact;
    PetscBool print_rollers; //
    PetscBool is_rotating;
 
    // PetscBool use_fieldsplit_for_bc;
 
    string contact_history;
    string move_history;
    string dirichlet_history;
    string mesh_file_str;
    string restart_file_str;
 
    boost::shared_ptr<BoundaryEle> update_roller;
    boost::shared_ptr<DomainEle> calc_reactions;
    boost::shared_ptr<DomainEle> contact_pipeline;
 
    moab::Core mb_post_debug;
    moab::Interface *moab_debug;
 
    int restart_step;
 
    Range skinEnts;
    std::map<EntityHandle, int> mapBlockTets;
    ProblemData() {
      order = 2;
      no_output = PETSC_FALSE;
      debug_info = PETSC_FALSE;
      output_freq = 1;
      atom_test_nb = 0;
      scale_params = PETSC_FALSE;
 
      is_large_strain = PETSC_FALSE;
      is_fieldsplit_bc_only = PETSC_FALSE;
 
      is_lambda_split = PETSC_FALSE;
 
      is_neohooke = PETSC_FALSE;
      with_plasticity = PETSC_FALSE;
 
      calculate_reactions = PETSC_FALSE;
      reaction_id = -1;
 
      with_contact = PETSC_FALSE;
      print_rollers = PETSC_FALSE; //
      is_rotating = PETSC_FALSE;
 
      // use_fieldsplit_for_bc = PETSC_TRUE;
 
      is_ho_geometry = PETSC_FALSE;
 
      save_restarts = PETSC_FALSE;
      is_restart = PETSC_FALSE;
 
      restart_step = -1;
 
      contact_history = move_history = dirichlet_history = "-load_history";
      moab_debug = &mb_post_debug;
    };
  };
 
  static ProblemData data;
 
  struct MMonitor : public FEMethod {
 
    MMonitor(
        SmartPetscObj<DM> &dm, MoFEM::Interface &m_field,
        boost::shared_ptr<OpContactTools::CommonData> common_data_ptr,
        boost::shared_ptr<PostProcVolumeOnRefinedMesh> post_proc_fe,
        std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> ux_scatter,
        std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uy_scatter,
        std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uz_scatter)
        : dM(dm), mField(m_field), commonDataPtr(common_data_ptr),
          postProcFe(post_proc_fe), uXScatter(ux_scatter),
          uYScatter(uy_scatter), uZScatter(uz_scatter) {
      volSideFe = boost::make_shared<MoFEM::DomainSideEle>(mField);
      volSideFe->getOpPtrVector().push_back(
          new OpVolumeSideCalculateTAU("TAU", commonDataPtr));
      volSideFe->getOpPtrVector().push_back(
          new OpVolumeSideCalculateEP("EP", commonDataPtr));
      postProcFeSkeleton =
          boost::make_shared<PostProcFaceOnRefinedMesh>(mField);
      postProcFeSkeleton->generateReferenceElementMesh();
      postProcFeSkeleton->getOpPtrVector().push_back(
          new OpCalculateJumpOnSkeleton("SKELETON_LAMBDA", commonDataPtr,
                                        volSideFe));
      postProcFeSkeleton->getOpPtrVector().push_back(
          new OpPostProcPlasticSkeleton(
              "SKELETON_LAMBDA", postProcFeSkeleton->postProcMesh,
              postProcFeSkeleton->mapGaussPts, commonDataPtr));
 
      postProcFeSkin = boost::make_shared<PostProcFaceOnRefinedMesh>(mField);
      CHKERR postProcFeSkin->generateReferenceElementMesh();
      CHKERR postProcFeSkin->addFieldValuesPostProc("U", "DISPLACEMENT");
      if (data.is_ho_geometry)
        postProcFeSkin->addFieldValuesPostProc("MESH_NODE_POSITIONS");
 
      if (data.with_plasticity) {
        CHKERR postProcFeSkin->addFieldValuesPostProcOnSkin("TAU");
        CHKERR postProcFeSkin->addFieldValuesPostProcOnSkin("EP");
      }
 
      if (data.with_contact) {
        CHKERR postProcFeSkin->addFieldValuesPostProcOnSkin("SIGMA");
      }
    };
 
    MoFEMErrorCode preProcess() {
      MoFEMFunctionBegin;
      CHKERR TSGetTimeStep(ts, &(commonDataPtr->timeStepSize));
      MoFEMFunctionReturn(0);
    }
    MoFEMErrorCode operator()() { return 0; }
 
    MoFEMErrorCode postProcess() {
      MoFEMFunctionBegin;
 
      auto make_vtk = [&]() {
        MoFEMFunctionBegin;
        if (data.no_output || ts_step % data.output_freq != 0)
          MoFEMFunctionReturnHot(0);
        CHKERR DMoFEMLoopFiniteElements(dM, "dFE", postProcFe);
        CHKERR postProcFe->writeFile("out_plastic_" +
                                     boost::lexical_cast<std::string>(ts_step) +
                                     ".h5m");
        MoFEMFunctionReturn(0);
      };
 
      auto save_skeleton = [&]() {
        MoFEMFunctionBegin;
 
        CHKERR DMoFEMLoopFiniteElements(dM, "sFE", postProcFeSkeleton);
        CHKERR postProcFeSkeleton->writeFile(
            "out_skeleton_" + boost::lexical_cast<std::string>(ts_step) +
            ".h5m");
        MoFEMFunctionReturn(0);
      };
 
      auto save_skin = [&]() {
        MoFEMFunctionBegin;
        if (data.no_output || ts_step % data.output_freq != 0)
          MoFEMFunctionReturnHot(0);
 
        CHKERR DMoFEMLoopFiniteElements(dM, "bFE", postProcFeSkin);
        CHKERR postProcFeSkin->writeFile(
            "out_skin_" + boost::lexical_cast<std::string>(ts_step) + ".h5m");
        MoFEMFunctionReturn(0);
      };
 
      double max_disp, min_disp;
 
      auto print_max_min = [&](auto &tuple, const std::string msg) {
        MoFEMFunctionBegin;
        CHKERR VecScatterBegin(std::get<1>(tuple), ts_u, std::get<0>(tuple),
                               INSERT_VALUES, SCATTER_FORWARD);
        CHKERR VecScatterEnd(std::get<1>(tuple), ts_u, std::get<0>(tuple),
                             INSERT_VALUES, SCATTER_FORWARD);
        CHKERR VecMax(std::get<0>(tuple), PETSC_NULL, &max_disp);
        CHKERR VecMin(std::get<0>(tuple), PETSC_NULL, &min_disp);
        PetscPrintf(PETSC_COMM_WORLD, "%s time %3.4e min %3.4e max %3.4e\n",
                    msg.c_str(), ts_t, min_disp, max_disp);
        MoFEMFunctionReturn(0);
      };
 
      double gauss_area;
 
      if (data.with_contact) {
        auto &l_reactions = commonDataPtr->bodyReactions;
        auto &lgauss_pts = commonDataPtr->stateArrayArea;
 
        int roller_nb = (*cache).rollerDataVec.size();
        l_reactions.resize(roller_nb * 3);
        std::fill(l_reactions.begin(), l_reactions.end(), 0);
 
        CHKERR DMoFEMLoopFiniteElements(dM, "bFE", data.update_roller);
 
        if (true) {
          std::vector<double> gauss_pts(2, 0);
          auto dm = mField.getInterface<Simple>()->getDM();
          CHKERR MPIU_Allreduce(lgauss_pts.data(), gauss_pts.data(), 2,
                                MPIU_REAL, MPIU_SUM,
                                PetscObjectComm((PetscObject)dm));
          MOFEM_LOG_C("WORLD", Sev::inform,
                      "\t TS %d Total contact area (ref): %g / %g", ts_step,
                      gauss_pts[0], gauss_pts[1]);
          gauss_area = gauss_pts[0];
          lgauss_pts[0] = lgauss_pts[1] = 0;
 
          vector<double> g_reactions(l_reactions.size() * 3, 0);
          CHKERR MPIU_Allreduce(l_reactions.data(), g_reactions.data(),
                                roller_nb * 3, MPIU_REAL, MPIU_SUM,
                                PetscObjectComm((PetscObject)dm));
          for (int dd = 0; dd != roller_nb; ++dd)
            MOFEM_LOG_C("WORLD", Sev::inform,
                        "\t Body %d reactions Time %g Force X: %3.4e Y: "
                        "%3.4e Z: %3.4e",
                        (*cache).rollerDataVec[dd].iD, ts_t,
                        g_reactions[3 * dd + 0], g_reactions[3 * dd + 1],
                        g_reactions[3 * dd + 2]);
        }
      }
 
      if (ts_step % data.output_freq == 0 && !data.no_output) {
        std::ostringstream ostrm, ostrm2;
        ostrm << "out_debug_" << ts_step << ".vtk";
        ostrm2 << "out_debug_" << ts_step << ".h5m";
        if (data.with_contact && data.debug_info) {
          if (mField.get_comm_size() == 1)
            CHKERR data.moab_debug->write_file(ostrm.str().c_str(), "VTK", "");
          else
            CHKERR data.moab_debug->write_file(ostrm2.str().c_str(), "MOAB",
                                               "PARALLEL=WRITE_PART");
 
          data.moab_debug->delete_mesh();
        }
      }
 
      if (data.calculate_reactions) {
        auto &vec = commonDataPtr->reactionsVec;
        CHKERR VecZeroEntries(vec);
        CHKERR VecGhostUpdateBegin(vec, INSERT_VALUES, SCATTER_FORWARD);
        CHKERR VecGhostUpdateEnd(vec, INSERT_VALUES, SCATTER_FORWARD);
        CHKERR DMoFEMLoopFiniteElements(dM, "dFE", data.calc_reactions);
        Simple *simple = mField.getInterface<Simple>();
        Range ents;
        CHKERR mField.get_moab().get_entities_by_dimension(0, 0, ents, true);
 
        ParallelComm *pcomm =
            ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
        CHKERR pcomm->reduce_tags(commonDataPtr->reactionTag, MPI_SUM, ents);
 
        CHKERR make_vtk();
        CHKERR save_skin();
 
        // clear reactions tag
        double def_value = 0.;
        CHKERR mField.get_moab().tag_clear_data(commonDataPtr->reactionTag,
                                                ents, &def_value);
 
        // CHKERR VecGhostUpdateBegin(vec, ADD_VALUES, SCATTER_REVERSE);
        // CHKERR VecGhostUpdateEnd(vec, ADD_VALUES, SCATTER_REVERSE);
        CHKERR VecAssemblyBegin(vec);
        CHKERR VecAssemblyEnd(vec);
        CHKERR VecGhostUpdateBegin(vec, ADD_VALUES, SCATTER_REVERSE);
        CHKERR VecGhostUpdateEnd(vec, ADD_VALUES, SCATTER_REVERSE);
        CHKERR VecGhostUpdateBegin(vec, INSERT_VALUES, SCATTER_FORWARD);
        CHKERR VecGhostUpdateEnd(vec, INSERT_VALUES, SCATTER_FORWARD);
 
        const double *react_ptr;
        CHKERR VecGetArrayRead(vec, &react_ptr);
 
        PetscPrintf(PETSC_COMM_WORLD,
                    "Reactions time %3.4e X: %3.4e Y: %3.4e Z: %3.4e \n", ts_t,
                    react_ptr[0], react_ptr[1], react_ptr[2]);
 
        CHKERR VecRestoreArrayRead(vec, &react_ptr);
      } else {
        CHKERR make_vtk();
        CHKERR save_skin();
      }
      // FIXME: this does not work properly
      if (data.debug_info && data.is_rotating && data.with_plasticity)
        CHKERR save_skeleton();
 
      if (data.with_contact && data.print_rollers &&
          ts_step % data.output_freq == 0 && !data.no_output) {
        auto &m_field = postProcFe->mField;
        data.moab_debug->delete_mesh();
        if (m_field.get_comm_rank() == 0) {
          std::ostringstream ost;
          for (auto &rol : (*cache).rollerDataVec) {
            EntityHandle vertex;
            VectorDouble roller_coords = rol.BodyDispScaled;
            roller_coords += rol.originCoords;
 
            // ugly workaroud to get correct first values
            if (!rol.positionDataParamName.empty() &&
                (ts_step == 0 || ts_step == data.restart_step)) {
              auto scale_method = boost::make_shared<TimeAccelerogram>(
                  rol.positionDataParamName);
              CHKERR MethodForForceScaling::applyScale(this, scale_method,
                                                       rol.BodyDispScaled);
              roller_coords = rol.BodyDispScaled;
            }
 
            CHKERR data.moab_debug->create_vertex(&roller_coords(0), vertex);
            CHKERR rol.saveBasicDataOnTag(*data.moab_debug, vertex);
          }
 
          ost << "out_roller_" << ts_step << ".vtk";
          CHKERR data.moab_debug->write_file(ost.str().c_str(), "VTK", "");
          data.moab_debug->delete_mesh();
        }
      }
 
      if (data.save_restarts && ts_step % data.output_freq == 0) {
 
        CHKERR VecGhostUpdateBegin(ts_u, INSERT_VALUES, SCATTER_FORWARD);
        CHKERR VecGhostUpdateEnd(ts_u, INSERT_VALUES, SCATTER_FORWARD);
 
        string rest_name =
            "restart_" + to_string(ts_step) + "_" + to_string(ts_t) + ".dat";
        PetscViewer viewer;
        CHKERR PetscViewerBinaryOpen(PETSC_COMM_WORLD, rest_name.c_str(),
                                     FILE_MODE_WRITE, &viewer);
        CHKERR PetscViewerBinarySetSkipInfo(viewer, PETSC_TRUE);
        VecView(ts_u, viewer);
        PetscViewerDestroy(&viewer);
      }
 
      CHKERR print_max_min(uXScatter, "Ux");
      CHKERR print_max_min(uYScatter, "Uy");
      CHKERR print_max_min(uZScatter, "Uz");
 
      switch (data.atom_test_nb) {
      case 1: {
        if (ts_step == 10) {
          const double *react_ptr;
          CHKERR VecGetArrayRead(commonDataPtr->reactionsVec, &react_ptr);
          if (fabs(react_ptr[0] + 0.12519621572) > 1e-6)
            SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
                    "atom test diverged!");
          CHKERR VecRestoreArrayRead(commonDataPtr->reactionsVec, &react_ptr);
        }
        break;
      }
      case 2: {
        if (ts_step == 5) {
          const double *react_ptr;
          CHKERR VecGetArrayRead(commonDataPtr->reactionsVec, &react_ptr);
          // hertz solution 0.25*(4/3)*(1/(1-0.3*0.3))*(10^0.5)*0.1^(3/2)
          if ((0.03663003663 - react_ptr[2]) / 0.03663003663 > 1e-2)
            SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
                    "atom test diverged!");
          CHKERR VecRestoreArrayRead(commonDataPtr->reactionsVec, &react_ptr);
        }
        break;
      }
      case 3: {
        if (ts_step == 50) {
          double min_disp;
          CHKERR VecMin(std::get<0>(uXScatter), PETSC_NULL, &min_disp);
          if (fabs(min_disp + 2.6949) > 1e-2)
            SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
                    "atom test diverged!");
        }
        break;
      }
      case 4: {
        if (ts_step == 5) {
          if (fabs(sqrt(gauss_area * 4.0 / M_PI) - 0.4574) / 0.4574 > 4e-2)
            SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
                    "atom test diverged!");
        }
        break;
      }
 
      default:
        break;
      }
 
      MoFEMFunctionReturn(0);
    }
 
  private:
    SmartPetscObj<DM> dM;
    MoFEM::Interface &mField;
    boost::shared_ptr<OpContactTools::CommonData> commonDataPtr;
    boost::shared_ptr<PostProcVolumeOnRefinedMesh> postProcFe;
 
    boost::shared_ptr<DomainSideEle> volSideFe;
    boost::shared_ptr<PostProcFaceOnRefinedMesh> postProcFeSkeleton;
    boost::shared_ptr<PostProcFaceOnRefinedMesh> postProcFeSkin;
 
    std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uXScatter;
    std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uYScatter;
    std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uZScatter;
  };
};
 
ContactPlasticity::ProblemData ContactPlasticity::data;
 
MoFEMErrorCode ContactPlasticity::runProblem() {
  MoFEMFunctionBegin;
  CHKERR setUP();
  CHKERR createCommonData();
  CHKERR bC();
  CHKERR OPs();
  CHKERR tsSolve();
  CHKERR postProcess();
  CHKERR checkResults();
  MoFEMFunctionReturn(0);
}
MoFEMErrorCode ContactPlasticity::getAnalysisParameters() {
  MoFEMFunctionBegin;
 
  CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "", "none");
  // Set approximation order
  CHKERR PetscOptionsInt("-order", "approximation order", "", data.order,
                         &data.order, PETSC_NULL);
  CHKERR PetscOptionsBool("-no_output", "save no output files", "",
                          data.no_output, &data.no_output, PETSC_NULL);
  CHKERR PetscOptionsBool("-debug_info", "print debug info", "",
                          data.debug_info, &data.debug_info, PETSC_NULL);
  CHKERR PetscOptionsBool("-scale_params",
                          "scale parameters (young modulus = 1)", "",
                          data.scale_params, &data.scale_params, PETSC_NULL);
  CHKERR PetscOptionsInt("-output_every",
                         "frequncy how often results are dumped on hard disk",
                         "", 1, &data.output_freq, PETSC_NULL);
  CHKERR PetscOptionsBool("-calculate_reactions", "calculate reactions", "",
                          data.calculate_reactions, &data.calculate_reactions,
                          PETSC_NULL);
  CHKERR PetscOptionsInt("-reaction_id", "meshset id for computing reactions",
                         "", data.reaction_id, &data.reaction_id, PETSC_NULL);
  CHKERR PetscOptionsInt("-atom_test", "number of atom test", "",
                         data.atom_test_nb, &data.atom_test_nb, PETSC_NULL);
  CHKERR PetscOptionsBool("-is_large_strain", "is large strains regime", "",
                          data.is_large_strain, &data.is_large_strain,
                          PETSC_NULL);
  CHKERR PetscOptionsBool(
      "-is_fieldsplit_bc_only", "use fieldsplit for boundary only", "",
      data.is_fieldsplit_bc_only, &data.is_fieldsplit_bc_only, PETSC_NULL);
  CHKERR PetscOptionsBool(
      "-is_lambda_split",
      "use axiator split approach, with reduced integration", "",
      data.is_lambda_split, &data.is_lambda_split, PETSC_NULL);
 
  CHKERR PetscOptionsBool("-is_neohooke", "is neohooke material", "",
                          data.is_neohooke, &data.is_neohooke, PETSC_NULL);
  // CHKERR PetscOptionsBool(
  //     "-use_fieldsplit_for_bc", "use fieldsplit for boundary conditions", "",
  //     data.use_fieldsplit_for_bc, &data.use_fieldsplit_for_bc, PETSC_NULL);
 
  CHKERR PetscOptionsBool(
      "-with_plasticity", "run calculations with plasticity", "",
      data.with_plasticity, &data.with_plasticity, PETSC_NULL);
 
  CHKERR PetscOptionsBool("-is_rotating", "is rotating frame used", "",
                          data.is_rotating, &data.is_rotating, PETSC_NULL);
  CHKERR PetscOptionsBool("-with_contact", "run calculations with contact", "",
                          data.with_contact, &data.with_contact, PETSC_NULL);
  CHKERR PetscOptionsBool("-print_rollers",
                          "output file with rollers positions", "",
                          data.print_rollers, &data.print_rollers, PETSC_NULL);
 
  char load_hist_file[255];
  PetscBool ctg_flag = PETSC_FALSE;
  CHKERR PetscOptionsString("-contact_history", "load_history.in", "",
                            "load_history.in", load_hist_file, 255, &ctg_flag);
  if (ctg_flag)
    data.contact_history = "-contact_history";
  CHKERR PetscOptionsString("-move_history", "load_history.in", "",
                            "load_history.in", load_hist_file, 255, &ctg_flag);
 
  if (ctg_flag)
    data.move_history = "-move_history";
 
  CHKERR PetscOptionsString("-dirichlet_history", "load_history.in", "",
                            "load_history.in", load_hist_file, 255, &ctg_flag);
  if (ctg_flag)
    data.dirichlet_history = "-dirichlet_history";
 
  char mesh_file_name[255];
  CHKERR PetscOptionsString("-file_name", "file name", "", "mesh.h5m",
                            mesh_file_name, 255, PETSC_NULL);
  data.mesh_file_str = string(mesh_file_name);
 
  char restart_file_name[255];
  CHKERR PetscOptionsString("-restart_file", "restart file name", "",
                            "restart.dat", restart_file_name, 255,
                            &data.is_restart);
  data.restart_file_str = string(restart_file_name);
 
  CHKERR PetscOptionsBool("-save_restarts",
                          "save restart files at each iteration", "",
                          data.save_restarts, &data.save_restarts, PETSC_NULL);
 
  ierr = PetscOptionsEnd();
  CHKERRQ(ierr);
 
  if (data.is_large_strain && data.with_plasticity && data.is_neohooke)
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Neohooke material is not supported with plasticity.");
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode ContactPlasticity::myMeshPartition(MoFEM::Interface &m_field) {
  MoFEMFunctionBegin;
 
  auto &moab = m_field.get_moab();
  Range tets;
  CHKERR moab.get_entities_by_dimension(0, 3, tets, false);
  Range faces;
  CHKERR moab.get_adjacencies(tets, 2, true, faces, moab::Interface::UNION);
  Range edges;
  CHKERR moab.get_adjacencies(tets, 1, true, edges, moab::Interface::UNION);
 
  CHKERR m_field.getInterface<ProblemsManager>()->partitionMesh(
      tets, 3, 0, m_field.get_comm_size());
 
  EntityHandle part_set;
  ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
  if (pcomm == NULL)
    pcomm = new ParallelComm(&moab, PETSC_COMM_WORLD);
  CHKERR moab.create_meshset(MESHSET_SET, part_set);
  Tag part_tag = pcomm->part_tag();
  Range proc_ents;
  Range all_proc_ents;
  Range off_proc_ents;
  Range tagged_sets;
  CHKERR moab.get_entities_by_type_and_tag(0, MBENTITYSET, &part_tag, NULL, 1,
                                           tagged_sets, moab::Interface::UNION);
 
  for (auto &mit : tagged_sets) {
    int part;
    CHKERR moab.tag_get_data(part_tag, &mit, 1, &part);
    if (part == m_field.get_comm_rank()) {
      CHKERR moab.get_entities_by_dimension(mit, 3, proc_ents, true);
      CHKERR moab.get_entities_by_handle(mit, all_proc_ents, true);
    } else
      CHKERR moab.get_entities_by_handle(mit, off_proc_ents, true);
  }
 
  Skinner skin(&moab);
  Range proc_ents_skin[4];
  proc_ents_skin[3] = proc_ents;
 
  Range all_tets, all_skin;
  CHKERR m_field.get_moab().get_entities_by_dimension(0, 3, all_tets, false);
  CHKERR skin.find_skin(0, all_tets, false, all_skin);
  CHKERR skin.find_skin(0, proc_ents, false, proc_ents_skin[2]);
 
  proc_ents_skin[2] = subtract(proc_ents_skin[2], all_skin);
  CHKERR moab.get_adjacencies(proc_ents_skin[2], 1, false, proc_ents_skin[1],
                              moab::Interface::UNION);
  CHKERR moab.get_connectivity(proc_ents_skin[1], proc_ents_skin[0], true);
  for (int dd = 0; dd != 3; dd++)
    CHKERR moab.add_entities(part_set, proc_ents_skin[dd]);
  CHKERR moab.add_entities(part_set, all_proc_ents);
 
  // not sure if necessary...
  {
    Range all_ents;
    CHKERR moab.get_entities_by_handle(0, all_ents);
    std::vector<unsigned char> pstat_tag(all_ents.size(), 0);
    CHKERR moab.tag_set_data(pcomm->pstatus_tag(), all_ents,
                             &*pstat_tag.begin());
  }
 
  Range ents_to_keep;
  CHKERR moab.get_entities_by_handle(part_set, ents_to_keep, false);
  off_proc_ents = subtract(off_proc_ents, ents_to_keep);
  Range meshsets;
  CHKERR moab.get_entities_by_type(0, MBENTITYSET, meshsets, false);
  for (auto m : meshsets) {
    CHKERR moab.remove_entities(m, off_proc_ents);
    // CHKERR moab.add_entities(part_set, &m, 1);
  }
  CHKERR moab.delete_entities(off_proc_ents);
 
  CHKERR pcomm->resolve_shared_ents(0, proc_ents, 3, -1, proc_ents_skin);
 
  BitRefLevel bit_level0;
  bit_level0.set(0);
  Simple *simple = m_field.getInterface<Simple>();
 
  Range ents;
  CHKERR m_field.get_moab().get_entities_by_dimension(part_set, 3, ents, true);
  CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevel(
      ents, simple->getBitRefLevel(), false);
  simple->getMeshSet() = part_set;
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode ContactPlasticity::loadMeshBlockSets() {
  MoFEMFunctionBegin;
  string block_name = "MAT_ELASTIC";
  // TODO: implement blocks for plasticity
  if (!data.with_plasticity)
    for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(mField, BLOCKSET, it)) {
      if (it->getName().compare(0, block_name.length(), block_name) == 0) {
        std::vector<double> block_data;
        CHKERR it->getAttributes(block_data);
        const int id = it->getMeshsetId();
        EntityHandle meshset = it->getMeshset();
        Range tets;
        // on a particular part, the meshset can be empty
        rval =
            mField.get_moab().get_entities_by_dimension(meshset, 3, tets, true);
        rval = MB_SUCCESS;
        for (auto ent : tets)
          data.mapBlockTets[ent] = id;
 
        mat_blocks[id] = BlockParamData();
        mat_blocks.at(id).young_modulus = block_data[0];
        mat_blocks.at(id).poisson = block_data[1];
        if (block_data.size() > 2)
          mat_blocks.at(id).density = block_data[2];
        mat_blocks.at(id).cn_cont = 1. / block_data[0];
 
        if (block_data.size() < 2)
          SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                  "provide an appropriate number of entries (2) parameters for "
                  "MAT_ELASTIC block");
      }
    }
 
  if (mat_blocks.size() < 2)
    mat_blocks[0] = BlockParamData();
 
  cache = &mat_blocks.begin()->second;
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode ContactPlasticity::setUP() {
  MoFEMFunctionBegin;
  Simple *simple = mField.getInterface<Simple>();
 
  // Select base
  // Note: For tets we have only H1 Ainsworth base, for Hex we have only H1
  // Demkowicz base. We need to implement Demkowicz H1 base on tet.
  enum bases { AINSWORTH, DEMKOWICZ, BERNSTEIN, LASBASETOPT };
  const char *list_bases[LASBASETOPT] = {"ainsworth", "demkowicz", "bernstein"};
  PetscInt choice_base_value = AINSWORTH;
  CHKERR PetscOptionsGetEList(PETSC_NULL, NULL, "-base", list_bases,
                              LASBASETOPT, &choice_base_value, PETSC_NULL);
 
  FieldApproximationBase base;
  switch (choice_base_value) {
  case AINSWORTH:
    base = AINSWORTH_LEGENDRE_BASE;
    MOFEM_LOG("WORLD", Sev::inform)
        << "Set AINSWORTH_LEGENDRE_BASE for displacements";
    break;
  case DEMKOWICZ:
    base = DEMKOWICZ_JACOBI_BASE;
    MOFEM_LOG("WORLD", Sev::inform)
        << "Set DEMKOWICZ_JACOBI_BASE for displacements";
    break;
  case BERNSTEIN:
    base = AINSWORTH_BERNSTEIN_BEZIER_BASE;
    MOFEM_LOG("WORLD", Sev::inform)
        << "Set AINSWORTH_BERNSTEIN_BEZIER_BASE for displacements";
    break;
  default:
    base = LASTBASE;
    break;
  }
 
  auto skin_ents = getEntsOnMeshSkin();
  data.skinEnts = skin_ents;
 
  Range tets, nodes;
  int nb_tets, nb_hexes;
  CHKERR mField.get_moab().get_entities_by_dimension(0, 3, tets);
  CHKERR mField.get_moab().get_connectivity(&*tets.begin(), 1, nodes);
  CHKERR mField.get_moab().get_number_entities_by_type(0, MBTET, nb_tets, true);
  CHKERR mField.get_moab().get_number_entities_by_type(0, MBHEX, nb_hexes,
                                                       true);
 
  if ((nb_tets && nodes.size() > 4) || (nb_hexes && nodes.size() > 8)) {
    data.is_ho_geometry = PETSC_TRUE;
    CHKERR simple->addDataField("MESH_NODE_POSITIONS", H1, base, 3);
    CHKERR simple->setFieldOrder("MESH_NODE_POSITIONS", 2);
    // auto all_edges_ptr = boost::make_shared<Range>();
    // CHKERR mField.get_moab().get_entities_by_type(0, MBEDGE, *all_edges_ptr,
    //                                               true);
    // CHKERR simple->setFieldOrder("MESH_NODE_POSITIONS", 2,
    // all_edges_ptr.get());
  }
  // Add field
  CHKERR simple->addDomainField("U", H1, base, 3);
  CHKERR simple->addBoundaryField("U", H1, base, 3);
  CHKERR simple->setFieldOrder("U", data.order);
 
  if (data.with_contact) {
    CHKERR simple->addDomainField("SIGMA", HDIV, DEMKOWICZ_JACOBI_BASE, 3);
    CHKERR simple->addBoundaryField("SIGMA", HDIV, DEMKOWICZ_JACOBI_BASE, 3);
    CHKERR simple->setFieldOrder("SIGMA", 0);
    CHKERR simple->setFieldOrder("SIGMA", data.order - 1, &skin_ents);
  }
  if (data.with_plasticity) {
    CHKERR simple->addDomainField("TAU", space_test, base, 1);
    CHKERR simple->setFieldOrder("TAU", std::max(0, data.order - 1));
    CHKERR simple->addDomainField("EP", space_test, base, 6);
    CHKERR simple->setFieldOrder("EP", std::max(0, data.order - 1));
  }
 
  if (data.is_rotating && data.with_plasticity && false) {
    CHKERR simple->addSkeletonField("SKELETON_LAMBDA", L2, base, 7);
    CHKERR simple->setFieldOrder("SKELETON_LAMBDA",
                                 std::max(0, data.order - 1));
  }
 
  CHKERR simple->defineFiniteElements();
 
  // Add Neumann forces elements
  CHKERR MetaNeumannForces::addNeumannBCElements(mField, "U");
  CHKERR MetaNodalForces::addElement(mField, "U");
  CHKERR MetaEdgeForces::addElement(mField, "U");
 
  simple->getOtherFiniteElements().push_back("FORCE_FE");
  simple->getOtherFiniteElements().push_back("PRESSURE_FE");
  if (data.is_rotating && data.with_plasticity && false)
    CHKERR simple->setSkeletonAdjacency();
 
  // PetscBool is_partitioned = PETSC_TRUE;
  CHKERR simple->defineProblem(PETSC_TRUE);
  CHKERR simple->buildFields();
 
  if (data.is_ho_geometry) {
    Projection10NodeCoordsOnField ent_method_material(mField,
                                                      "MESH_NODE_POSITIONS");
    CHKERR mField.loop_dofs("MESH_NODE_POSITIONS", ent_method_material);
  }
 
  CHKERR simple->buildFiniteElements();
  if (data.is_rotating && data.with_plasticity && false) {
    CHKERR mField.remove_ents_from_finite_element("sFE", skin_ents);
    CHKERR mField.build_finite_elements("sFE");
  }
  CHKERR simple->buildProblem();
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode ContactPlasticity::createCommonData() {
  MoFEMFunctionBegin;
 
  auto get_material_stiffness = [&]() {
    MoFEMFunctionBegin;
 
    constexpr auto t_kd = Kronecker_Delta<int>();
    constexpr auto t_kds = Kronecker_Delta_symmetric<int>();
 
    // auto &t_D = commonDataPtr->tD;
    (*cache).Is(i, j, k, l) =
        (0.5 * t_kd(i, k) * t_kd(j, l)) + (0.5 * t_kd(i, l) * t_kd(j, k));
 
    // Tensor4<double, 3, 3, 3, 3> II;
    // II(i, j, k, l) = t_kd(i, k) * t_kd(j, l);
 
    auto Is_ddg = tensor4_to_ddg((*cache).Is);
 
    for (auto &mit : mat_blocks) {
 
      cache = &mit.second;
      const double G = (*cache).young_modulus / (2. * (1. + (*cache).poisson));
      const double K =
          (*cache).young_modulus / (3. * (1. - 2. * (*cache).poisson));
 
      // for plane stress
      double A = 6. * G / (3. * K + 4. * G);
      // is_plane_strain and 3D
      A = 1.;
 
      auto t_D_tmp = getFTensor4DdgFromMat<3, 3, 0>(*commonDataPtr->mtD);
      auto t_D_axiator =
          getFTensor4DdgFromMat<3, 3, 0>(*commonDataPtr->mtD_Axiator);
      auto t_D_deviator =
          getFTensor4DdgFromMat<3, 3, 0>(*commonDataPtr->mtD_Deviator);
 
      t_D_axiator(i, j, k, l) =
          A * (K - (2. / 3.) * G) * t_kds(i, j) * t_kds(k, l);
      t_D_deviator(i, j, k, l) = 2 * G * ((t_kds(i, k) ^ t_kds(j, l)) / 4.);
 
      t_D_tmp(i, j, k, l) = t_D_deviator(i, j, k, l) + t_D_axiator(i, j, k, l);
 
      (*cache).mtD = *commonDataPtr->mtD;
      (*cache).scale_constraint = 1.;
    }
 
    // commonDataPtr->isDualBase = data.is_dual_base;
 
    MoFEMFunctionReturn(0);
  };
 
  commonDataPtr = boost::make_shared<OpContactTools::CommonData>();
 
  commonDataPtr->mtD = boost::make_shared<MatrixDouble>();
  commonDataPtr->mtD->resize(36, 1);
  commonDataPtr->mtD_Axiator = boost::make_shared<MatrixDouble>();
  commonDataPtr->mtD_Axiator->resize(36, 1);
  commonDataPtr->mtD_Deviator = boost::make_shared<MatrixDouble>();
  commonDataPtr->mtD_Deviator->resize(36, 1);
 
  commonDataPtr->mGradPtr = boost::make_shared<MatrixDouble>();
  commonDataPtr->mStrainPtr = boost::make_shared<MatrixDouble>();
  commonDataPtr->mStressPtr = boost::make_shared<MatrixDouble>();
  commonDataPtr->mPiolaStressPtr = boost::make_shared<MatrixDouble>();
  commonDataPtr->mPiolaStressPtr->resize(9, 1, false);
 
  commonDataPtr->mDeformationGradient = boost::make_shared<MatrixDouble>();
  commonDataPtr->matC = boost::make_shared<MatrixDouble>();
  commonDataPtr->matEigenVal = boost::make_shared<MatrixDouble>();
  commonDataPtr->matEigenVector = boost::make_shared<MatrixDouble>();
  commonDataPtr->detF = boost::make_shared<VectorDouble>();
  commonDataPtr->materialTangent = boost::make_shared<MatrixDouble>();
  commonDataPtr->dE_dF_mat = boost::make_shared<MatrixDouble>();
  commonDataPtr->dE_dF_mat->resize(81, 1, false);
 
  commonDataPtr->contactStressPtr = boost::make_shared<MatrixDouble>();
  commonDataPtr->contactStressDivergencePtr =
      boost::make_shared<MatrixDouble>();
  commonDataPtr->contactNormalPtr = boost::make_shared<MatrixDouble>();
  commonDataPtr->contactGapPtr = boost::make_shared<VectorDouble>();
  commonDataPtr->contactGapDiffPtr = boost::make_shared<MatrixDouble>();
  commonDataPtr->contactNormalDiffPtr = boost::make_shared<MatrixDouble>();
 
  commonDataPtr->contactTractionPtr = boost::make_shared<MatrixDouble>();
  commonDataPtr->mDispPtr = boost::make_shared<MatrixDouble>();
 
  commonDataPtr->plasticSurfacePtr = boost::make_shared<VectorDouble>();
  commonDataPtr->plasticFlowPtr = boost::make_shared<MatrixDouble>();
  commonDataPtr->plasticTauPtr = boost::make_shared<VectorDouble>();
  commonDataPtr->plasticStrainPtr = boost::make_shared<MatrixDouble>();
 
  commonDataPtr->plasticTauDotPtr = boost::make_shared<VectorDouble>();
  commonDataPtr->plasticStrainDotPtr = boost::make_shared<MatrixDouble>();
 
  commonDataPtr->plasticTauJumpPtr = boost::make_shared<VectorDouble>();
  commonDataPtr->plasticStrainJumpPtr = boost::make_shared<MatrixDouble>();
  commonDataPtr->guidingVelocityPtr = boost::make_shared<MatrixDouble>();
  commonDataPtr->guidingVelocityPtr->resize(3, 1, false);
 
  commonDataPtr->plasticGradTauPtr = boost::make_shared<MatrixDouble>();
  commonDataPtr->plasticGradStrainPtr = boost::make_shared<MatrixDouble>();
 
  jAc.resize(3, 3, false);
  invJac.resize(3, 3, false);
 
  if (data.with_plasticity)
    CHKERR commonDataPtr->calculateDiffDeviator();
 
  CHKERR get_material_stiffness();
  if (data.is_neohooke)
    commonDataPtr->isNeoHook = true;
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode ContactPlasticity::bC() {
  MoFEMFunctionBegin;
  auto bc_mng = mField.getInterface<BcManager>();
  auto prb_mng = mField.getInterface<ProblemsManager>();
  auto simple = mField.getInterface<Simple>();
 
  auto fix_disp = [&](const std::string blockset_name) {
    Range fix_ents;
    for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(mField, BLOCKSET, it)) {
      if (it->getName().compare(0, blockset_name.length(), blockset_name) ==
          0) {
        CHKERR mField.get_moab().get_entities_by_handle(it->meshset, fix_ents,
                                                        true);
      }
    }
    return fix_ents;
  };
 
  auto remove_ents = [&](string prob, const Range &&ents, const int low,
                         const int high) {
    auto prb_mng = mField.getInterface<ProblemsManager>();
    auto simple = mField.getInterface<Simple>();
    MoFEMFunctionBegin;
    Range verts;
    CHKERR mField.get_moab().get_connectivity(ents, verts, true);
    CHKERR mField.get_moab().get_adjacencies(ents, 1, false, verts,
                                             moab::Interface::UNION);
    verts.merge(ents);
    CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(verts);
    CHKERR prb_mng->removeDofsOnEntities(prob, "U", verts, low, high);
    MoFEMFunctionReturn(0);
  };
 
  if (data.with_contact) {
    Range boundary_ents;
    boundary_ents.merge(fix_disp("FIX_X"));
    boundary_ents.merge(fix_disp("FIX_Y"));
    boundary_ents.merge(fix_disp("FIX_Z"));
    boundary_ents.merge(fix_disp("FIX_ALL"));
    boundary_ents.merge(fix_disp("NO_CONTACT"));
    CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
        boundary_ents);
    // CHKERR mField.getInterface<ProblemsManager>()->removeDofsOnEntities(
    //     mField.getInterface<Simple>()->getProblemName(), "SIGMA", boundary_ents,
    //     0, 2);
  }
 
  CHKERR remove_ents(simple->getProblemName(), fix_disp("FIX_X"), 0, 0);
  CHKERR remove_ents(simple->getProblemName(), fix_disp("FIX_Y"), 1, 1);
  CHKERR remove_ents(simple->getProblemName(), fix_disp("FIX_Z"), 2, 2);
  CHKERR remove_ents(simple->getProblemName(), fix_disp("FIX_ALL"), 0, 2);
 
  std::vector<DataFromBc> bcData;
  dirichlet_bc_ptr = boost::make_shared<DirichletDisplacementRemoveDofsBc>(
      mField, "U", simple->getProblemName(), "DISPLACEMENT", true);
  dirichlet_bc_ptr->iNitialize(); // TODO: check whether to use preProcess
                                  // better
 
  auto dm = simple->getDM();
  if (data.with_contact && data.with_plasticity) {
    CHKERR set_contact_dm(mField, dm, dmContact);
    CHKERR set_plastic_ep_dm(mField, dmContact, dmEP);
    CHKERR set_plastic_tau_dm(mField, dmEP, dmTAU);
  } else if (data.with_plasticity) {
    CHKERR set_plastic_ep_dm(mField, dm, dmEP);
    CHKERR set_plastic_tau_dm(mField, dmEP, dmTAU);
  }
 
  // forces and pressures on surface
  CHKERR MetaNeumannForces::setMomentumFluxOperators(mField, neumann_forces,
                                                     PETSC_NULL, "U");
  // nodal forces
  CHKERR MetaNodalForces::setOperators(mField, nodal_forces, PETSC_NULL, "U");
  // edge forces
  CHKERR MetaEdgeForces::setOperators(mField, edge_forces, PETSC_NULL, "U");
 
  // FIXME: add higher order operators !!
  // TODO: add higher order operators (check nonlinear_dynamics)
  for (auto mit = neumann_forces.begin(); mit != neumann_forces.end(); mit++) {
    if (data.is_ho_geometry)
      CHKERR addHOOpsFace3D("MESH_NODE_POSITIONS", mit->second->getLoopFe(),
                            false, false);
    mit->second->methodsOp.push_back(
        new TimeForceScale("-load_history", false));
    mit->second->methodsOp.push_back(new LoadScale((*cache).young_modulus_inv));
  }
  for (auto mit = nodal_forces.begin(); mit != nodal_forces.end(); mit++) {
    if (data.is_ho_geometry)
      CHKERR addHOOpsFace3D("MESH_NODE_POSITIONS", mit->second->getLoopFe(),
                            false, false);
    mit->second->methodsOp.push_back(
        new TimeForceScale("-load_history", false));
    mit->second->methodsOp.push_back(new LoadScale((*cache).young_modulus_inv));
  }
  for (auto mit = edge_forces.begin(); mit != edge_forces.end(); mit++) {
    if (data.is_ho_geometry)
      CHKERR addHOOpsFace3D("MESH_NODE_POSITIONS", mit->second->getLoopFe(),
                            false, false);
    mit->second->methodsOp.push_back(
        new TimeForceScale("-load_history", false));
    mit->second->methodsOp.push_back(new LoadScale((*cache).young_modulus_inv));
  }
 
  // dirichlet_bc_ptr->dIag = 1;
  dirichlet_bc_ptr->methodsOp.push_back(
      new TimeForceScale(data.dirichlet_history, false));
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode ContactPlasticity::OPs() {
  MoFEMFunctionBegin;
  PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
  Simple *simple = mField.getInterface<Simple>();
  auto problem_manager = mField.getInterface<ProblemsManager>();
  const Field_multiIndex *fields_ptr;
  CHKERR mField.get_fields(&fields_ptr);
 
  auto get_boundary_markers = [&](string prob, Range bound_ents, int lo,
                                  int hi) {
    auto marker_ptr = boost::make_shared<std::vector<unsigned char>>();
 
    CHKERR problem_manager->modifyMarkDofs(
        prob, ROW, "U", lo, hi, ProblemsManager::MarkOP::OR, 1, *marker_ptr);
 
    CHKERR problem_manager->markDofs(prob, ROW, ProblemsManager::AND,
                                     bound_ents, *marker_ptr);
 
    return marker_ptr;
  };
 
  typedef boost::shared_ptr<std::vector<unsigned char>> MarkerPtr;
 
  std::vector<MarkerPtr> markers_vec;
  boundaryMarker = boost::make_shared<std::vector<char unsigned>>();
  allBoundaryEnts = boost::make_shared<Range>();
 
  auto get_adj = [&](Range &ents) {
    Range verts;
    CHKERR mField.get_moab().get_connectivity(ents, verts, true);
    CHKERR mField.get_moab().get_adjacencies(ents, 1, false, verts,
                                             moab::Interface::UNION);
    verts.merge(ents);
 
    CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(verts);
    return verts;
  };
 
  auto add_disp_bc_ops = [&]() {
    MoFEMFunctionBeginHot;
 
    if (data.with_plasticity)
      pipeline_mng->getOpBoundaryLhsPipeline().push_back(
          new OpSchurBeginBoundary("U", commonDataPtr));
 
    std::map<char, size_t> xyz_map = {{'X', 0}, {'Y', 1}, {'Z', 2}};
    // if(false)
    for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(mField, BLOCKSET, it)) {
      size_t ii = 0;
      for (string move_s : {"MOVE_X", "MOVE_Y", "MOVE_Z"}) {
        if (it->getName().compare(0, move_s.length(), move_s) == 0) {
          size_t idx =
              xyz_map.at(it->getName()[(it->getName().find("MOVE_") + 5)]);
          Range ents;
          CHKERR mField.get_moab().get_entities_by_handle(it->meshset, ents,
                                                          true);
          std::vector<double> disp_vec;
          it->getAttributes(disp_vec);
          if (disp_vec.empty())
            SETERRQ(
                PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                "provide an appropriate number of params (1) for MOVE block");
 
          VectorDouble disp({0, 0, 0});
          disp(ii) = disp_vec[0];
          auto verts = get_adj(ents);
          allBoundaryEnts->merge(verts);
          auto markers =
              get_boundary_markers(simple->getProblemName(), verts, ii, ii);
          markers_vec.push_back(markers);
 
          bcData.push_back(new DataFromMove(disp, verts));
          // rhs
          pipeline_mng->getOpBoundaryRhsPipeline().push_back(
              new OpSetBc("U", false, markers));
          pipeline_mng->getOpBoundaryRhsPipeline().push_back(
              new OpEssentialRHS_U("U", commonDataPtr,
                                   bcData.back().scaledValues,
                                   bcData.back().ents));
          pipeline_mng->getOpBoundaryRhsPipeline().push_back(
              new OpUnSetBc("U"));
          // lhs
          pipeline_mng->getOpBoundaryLhsPipeline().push_back(
              new OpSetBc("U", false, markers));
          pipeline_mng->getOpBoundaryLhsPipeline().push_back(
              new OpSpringLhsNoFS(
                  "U", "U", [](double, double, double) { return 1.; },
                  boost::make_shared<Range>(bcData.back().ents)));
          pipeline_mng->getOpBoundaryLhsPipeline().push_back(
              new OpUnSetBc("U"));
        }
        ++ii;
      }
    }
 
    auto merge_boundary_markers = [&](std::vector<MarkerPtr> &mark) {
      auto all_markers = boost::make_shared<std::vector<char unsigned>>();
      for (auto &mit : mark) {
        all_markers->resize(mit->size(), 0);
        for (int i = 0; i != mit->size(); ++i)
          (*all_markers)[i] |= (*mit)[i];
      }
      return all_markers;
    };
 
    boundaryMarker = merge_boundary_markers(markers_vec);
 
    auto bc_mng = mField.getInterface<BcManager>();
    CHKERR bc_mng->pushMarkDOFsOnEntities(simple->getProblemName(), "ROTATE_X",
                                          "U", 1, 2);
    // TODO:  implement this
    // CHKERR bc_mng->pushMarkDOFsOnEntities(simple->getProblemName(),
    // "ROTATE_Y",
    //                                       "U", 0, 0);
    CHKERR bc_mng->pushMarkDOFsOnEntities(simple->getProblemName(), "ROTATE_Z",
                                          "U", 0, 1);
    CHKERR bc_mng->pushMarkDOFsOnEntities(simple->getProblemName(),
                                          "ROTATE_ALL", "U", 0, 3);
 
    auto &bc_map = bc_mng->getBcMapByBlockName();
    if (bc_map.size()) {
      for (auto b : bc_map) {
        if (std::regex_match(b.first, std::regex("(.*)_ROTATE_(.*)"))) {
          boundaryMarker->resize(b.second->bcMarkers.size(), 0);
          for (int i = 0; i != b.second->bcMarkers.size(); ++i) {
            // CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
            //     b.second->bcEdges);
            allBoundaryEnts->merge(b.second->bcEnts);
            (*boundaryMarker)[i] |= b.second->bcMarkers[i];
          }
        }
      }
    }
 
    if (boundaryMarker->empty() && bc_map.size() == 0)
      boundaryMarker =
          get_boundary_markers(simple->getProblemName(), Range(), -1, -1);
 
    for (auto &bc : mField.getInterface<BcManager>()->getBcMapByBlockName()) {
      if (std::regex_match(bc.first, std::regex("(.*)_ROTATE_(.*)"))) {
        MOFEM_LOG("WORLD", Sev::inform)
            << "Set boundary (rotation) residual for " << bc.first;
        auto angles = boost::make_shared<VectorDouble>(3);
        auto c_coords = boost::make_shared<VectorDouble>(3);
        if (bc.second->bcAttributes.size() < 6)
          SETERRQ1(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
                   "Wrong size of boundary attributes vector. Set the correct "
                   "block size attributes (3) angles and (3) coordinates for "
                   "center of rotation. Size of attributes %d",
                   bc.second->bcAttributes.size());
        std::copy(&bc.second->bcAttributes[0], &bc.second->bcAttributes[3],
                  angles->data().begin());
        std::copy(&bc.second->bcAttributes[3], &bc.second->bcAttributes[6],
                  c_coords->data().begin());
        bcData.push_back(new DataFromMove(*angles, *bc.second->getBcEntsPtr()));
        // rhs
        pipeline_mng->getOpBoundaryRhsPipeline().push_back(
            new OpSetBc("U", false, bc.second->getBcMarkersPtr()));
        pipeline_mng->getOpBoundaryRhsPipeline().push_back(
            new OpRotate("U", bcData.back().scaledValues, c_coords,
                         bc.second->getBcEntsPtr()));
        pipeline_mng->getOpBoundaryRhsPipeline().push_back(new OpSpringRhs(
            "U", commonDataPtr->mDispPtr,
            [](double, double, double) { return 1.; },
            bc.second->getBcEntsPtr()));
        pipeline_mng->getOpBoundaryRhsPipeline().push_back(new OpUnSetBc("U"));
        // lhs
        pipeline_mng->getOpBoundaryLhsPipeline().push_back(
            new OpSetBc("U", false, bc.second->getBcMarkersPtr()));
        pipeline_mng->getOpBoundaryLhsPipeline().push_back(new OpSpringLhsNoFS(
            "U", "U", [](double, double, double) { return 1.; },
            bc.second->getBcEntsPtr()));
        pipeline_mng->getOpBoundaryLhsPipeline().push_back(new OpUnSetBc("U"));
      }
    }
 
    MoFEMFunctionReturnHot(0);
  };
 
  auto add_ho_ops = [&](auto &pipeline) {
    MoFEMFunctionBeginHot;
    auto jac_ptr = boost::make_shared<MatrixDouble>();
    auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
    auto det_ptr = boost::make_shared<VectorDouble>();
    pipeline.push_back(new OpCalculateVectorFieldGradient<3, 3>(
        "MESH_NODE_POSITIONS", jac_ptr));
    // pipeline.push_back(new OpCalculateHOJacVolume(jac_ptr));
    pipeline.push_back(new OpInvertMatrix<3>(jac_ptr, det_ptr, inv_jac_ptr));
    pipeline.push_back(
        new OpSetHOContravariantPiolaTransform(HDIV, det_ptr, jac_ptr));
    pipeline.push_back(new OpSetHOWeights(det_ptr));
    pipeline.push_back(new OpSetHOInvJacVectorBase(HDIV, inv_jac_ptr));
    pipeline.push_back(new OpSetHOInvJacToScalarBases(H1, inv_jac_ptr));
    pipeline.push_back(new OpSetHOWeights(det_ptr));
    // pipeline.push_back(new OpMakeHighOrderGeometryWeightsOnVolume());
    // pipeline.push_back(new OpSetHOGeometryCoordsOnVolume());
    MoFEMFunctionReturnHot(0);
  };
 
  auto add_domain_base_ops = [&](auto &pipeline) {
    MoFEMFunctionBeginHot;
    pipeline.push_back(
        new OpSetMaterialBlock("U", commonDataPtr, data.mapBlockTets));
    if (data.is_ho_geometry)
      CHKERR add_ho_ops(pipeline);
 
    pipeline.push_back(
        new OpCalculateVectorFieldValues<3>("U", commonDataPtr->mDispPtr));
    pipeline.push_back(
        new OpCalculateVectorFieldGradient<3, 3>("U", commonDataPtr->mGradPtr));
 
    if (data.with_plasticity) {
      pipeline.push_back(new OpCalculateTensor2SymmetricFieldValues<3>(
          "EP", commonDataPtr->plasticStrainPtr));
      pipeline.push_back(new OpCalculateScalarFieldValues(
          "TAU", commonDataPtr->plasticTauPtr));
      pipeline.push_back(new OpCalculateScalarFieldValuesDot(
          "TAU", commonDataPtr->plasticTauDotPtr));
      pipeline.push_back(new OpCalculateTensor2SymmetricFieldValuesDot<3>(
          "EP", commonDataPtr->plasticStrainDotPtr));
 
      if (data.is_rotating) {
        pipeline.push_back(new OpCalculateTensor2SymmetricFieldGradient<3, 3>(
            "EP", commonDataPtr->plasticGradStrainPtr));
        pipeline.push_back(new OpCalculateScalarFieldGradient<3>(
            "TAU", commonDataPtr->plasticGradTauPtr));
        pipeline.push_back(
            new OpCalculatePlasticConvRotatingFrame("TAU", commonDataPtr));
      }
    }
 
    if (data.is_large_strain) {
      pipeline.push_back(new OpLogStrain("U", commonDataPtr));
    } else
      pipeline.push_back(new OpStrain("U", commonDataPtr));
 
    MoFEMFunctionReturnHot(0);
  };
 
  auto add_domain_stress_ops = [&](auto &pipeline, auto mat_D_ptr) {
    MoFEMFunctionBeginHot;
 
    if (data.with_plasticity) {
      pipeline.push_back(new OpPlasticStress("U", commonDataPtr, mat_D_ptr));
 
      pipeline.push_back(
          new OpCalculatePlasticSurfaceAndFlow("U", commonDataPtr));
 
    } else
      pipeline.push_back(new OpStress("U", commonDataPtr, mat_D_ptr));
 
    MoFEMFunctionReturnHot(0);
  };
 
  auto add_skeleton_base_ops = [&](auto &pipeline) {
    volSideFe = boost::make_shared<MoFEM::DomainSideEle>(mField);
 
    volSideFe->getOpPtrVector().push_back(
        new OpVolumeSideCalculateTAU("TAU", commonDataPtr));
    volSideFe->getOpPtrVector().push_back(
        new OpVolumeSideCalculateEP("EP", commonDataPtr));
    // pipeline.push_back(new OpCalculateJumpOnSkeleton("SKELETON_LAMBDA",
    //                                                  commonDataPtr,
    //                                                  volSideFe));
  };
 
  auto add_domain_ops_lhs = [&](auto &pipeline, auto m_D_ptr) {
    MoFEMFunctionBeginHot;
 
    pipeline.push_back(new OpSetBc("U", true, boundaryMarker));
    if (data.with_plasticity)
      pipeline.push_back(new OpSchurBegin("U", commonDataPtr));
 
    if (data.is_large_strain) {
      if (data.is_neohooke)
        pipeline.push_back(
            new OpCalculateNeoHookeStressTangent("U", commonDataPtr));
      else
        pipeline.push_back(
            new OpCalculateLogStressTangent<true>("U", commonDataPtr, m_D_ptr));
      pipeline.push_back(
          new OpLogStrainMatrixLhs("U", "U", commonDataPtr->materialTangent));
    } else
      pipeline.push_back(new OpStiffnessMatrixLhs("U", "U", m_D_ptr));
    if (data.with_plasticity) {
      if (data.is_large_strain)
        pipeline.push_back(new OpCalculatePlasticInternalForceLhs_dEP<true>(
            "U", "EP", commonDataPtr, m_D_ptr));
      else
        pipeline.push_back(new OpCalculatePlasticInternalForceLhs_dEP<false>(
            "U", "EP", commonDataPtr, m_D_ptr));
      // FIXME: dummy one for testing
      // TODO: this can be deleted
      // pipeline.push_back(new OpCalculatePlasticInternalForceLhs_dTAU(
      //     "U", "TAU", commonDataPtr));
    }
 
    if (data.is_rotating)
      pipeline.push_back(new OpRotatingFrameLhs("U", "U", commonDataPtr));
 
    pipeline.push_back(new OpUnSetBc("U"));
 
    MoFEMFunctionReturnHot(0);
  };
 
  auto add_domain_ops_rhs = [&](auto &pipeline, auto m_D_ptr) {
    MoFEMFunctionBeginHot;
 
    auto get_body_force = [this](const double, const double, const double) {
      auto *pipeline_mng = mField.getInterface<PipelineManager>();
      auto fe_domain_rhs = pipeline_mng->getDomainRhsFE();
      Tensor1<double, 3> t_source;
      for (int i = 0; i != 3; i++)
        t_source(i) = (*cache).gravity[i] * (*cache).density;
      const auto time = fe_domain_rhs->ts_t;
      t_source(i) *= time;
      return t_source;
    };
    auto get_centrifugal_force = [this](const double x, const double y,
                                        const double z) {
      Tensor1<double, 3> t_coords(x, y, z);
      Tensor1<double, 3> t_source;
      t_source(i) = (*cache).density * (*cache).Omega(i, k) *
                    (*cache).Omega(k, j) * t_coords(j);
      // auto *pipeline_mng = mField.getInterface<PipelineManager>();
      // auto fe_domain_rhs = pipeline_mng->getDomainRhsFE();
      // const auto time = fe_domain_rhs->ts_t;
      // t_source(i) *= time;
      return t_source;
    };
 
    pipeline.push_back(new OpSetBc("U", true, boundaryMarker));
 
    if (data.is_rotating) {
      pipeline.push_back(new OpRotatingFrameRhs("U", commonDataPtr));
      // pipeline.push_back(new OpCentrifugalForce2("U",
      // get_centrifugal_force));
    }
 
    pipeline.push_back(new OpBodyForce("U", get_body_force));
    if (data.is_large_strain) {
      if (data.is_neohooke)
        pipeline.push_back(
            new OpCalculateNeoHookeStressTangent("U", commonDataPtr));
      else
        pipeline.push_back(new OpCalculateLogStressTangent<false>(
            "U", commonDataPtr, m_D_ptr));
      pipeline.push_back(
          new OpInternalForceRhs<true, false>("U", commonDataPtr));
    } else
      pipeline.push_back(
          new OpInternalForceRhs<false, false>("U", commonDataPtr));
 
    pipeline.push_back(new OpUnSetBc("U"));
    MoFEMFunctionReturnHot(0);
  };
 
  auto add_domain_ops_rhs_constraints = [&](auto &pipeline) {
    MoFEMFunctionBeginHot;
    if (data.with_plasticity) {
      pipeline.push_back(new OpCalculatePlasticFlowRhs("EP", commonDataPtr));
      pipeline.push_back(new OpCalculateConstraintRhs("TAU", commonDataPtr));
      if (data.debug_info)
        pipeline.push_back(new OpGetGaussPtsPlasticState("U", commonDataPtr));
    }
    MoFEMFunctionReturnHot(0);
  };
 
  auto add_domain_ops_lhs_constraints = [&](auto &pipeline, auto m_D_ptr) {
    MoFEMFunctionBeginHot;
 
    if (data.with_plasticity) {
 
      pipeline.push_back(new OpSetBc("U", true, boundaryMarker));
 
      if (data.is_large_strain) {
        pipeline.push_back(new OpCalculateLogStressTangent<false>(
            "U", commonDataPtr, m_D_ptr));
        pipeline.push_back(
            new OpCalculatePlasticFlowLhs_dU<true>("EP", "U", commonDataPtr));
        pipeline.push_back(
            new OpCalculateConstraintLhs_dU<true>("TAU", "U", commonDataPtr));
      } else {
        pipeline.push_back(
            new OpCalculatePlasticFlowLhs_dU<false>("EP", "U", commonDataPtr));
        pipeline.push_back(
            new OpCalculateConstraintLhs_dU<false>("TAU", "U", commonDataPtr));
      }
 
      pipeline.push_back(
          new OpCalculatePlasticFlowLhs_dEP("EP", "EP", commonDataPtr));
      pipeline.push_back(
          new OpCalculatePlasticFlowLhs_dTAU("EP", "TAU", commonDataPtr));
      pipeline.push_back(
          new OpCalculateConstraintLhs_dEP("TAU", "EP", commonDataPtr));
      pipeline.push_back(
          new OpCalculateConstraintLhs_dTAU("TAU", "TAU", commonDataPtr));
 
      if (data.with_plasticity) {
        pipeline.push_back(new OpSchurPreconditionMassTau(
            "TAU", "TAU", [](double, double, double) { return 1.; },
            commonDataPtr));
        pipeline.push_back(new OpSchurEnd("U", commonDataPtr));
      }
      pipeline.push_back(new OpUnSetBc("U"));
    }
 
    MoFEMFunctionReturnHot(0);
  };
 
  auto add_domain_ops_contact_lhs = [&](auto &pipeline) {
    if (data.with_contact) {
      pipeline.push_back(new OpSetBc("U", true, boundaryMarker));
 
      pipeline.push_back(new OpMixDivULhs(
          "SIGMA", "U", []() { return 1.; }, true));
      pipeline.push_back(new OpLambdaGraULhs(
          "SIGMA", "U", []() { return 1.; }, true));
 
      pipeline.push_back(new OpSchurPreconditionMassContactVol(
          "SIGMA", "SIGMA", [](double, double, double) { return 1.; },
          commonDataPtr));
 
      pipeline.push_back(new OpUnSetBc("U"));
    }
  };
 
  auto add_domain_ops_contact_rhs = [&](auto &pipeline) {
    if (data.with_contact) {
      pipeline.push_back(new OpCalculateHVecTensorField<3, 3>(
          "SIGMA", commonDataPtr->contactStressPtr));
      pipeline.push_back(new OpCalculateHVecTensorDivergence<3, 3>(
          "SIGMA", commonDataPtr->contactStressDivergencePtr));
      pipeline.push_back(
          new OpMixDivURhs("SIGMA", commonDataPtr->mDispPtr,
                           [](double, double, double) { return 1; }));
      pipeline.push_back(
          new OpMiXLambdaGradURhs("SIGMA", commonDataPtr->mGradPtr));
 
      pipeline.push_back(new OpSetBc("U", true, boundaryMarker));
      pipeline.push_back(new OpMixUTimesDivLambdaRhs(
          "U", commonDataPtr->contactStressDivergencePtr));
      pipeline.push_back(
          new OpMixUTimesLambdaRhs("U", commonDataPtr->contactStressPtr));
      pipeline.push_back(new OpUnSetBc("U"));
    }
  };
 
  auto add_boundary_base_ops = [&](auto &pipeline) {
    MoFEMFunctionBeginHot;
    pipeline.push_back(
        new OpSetMaterialBlockBoundary("U", commonDataPtr, data.mapBlockTets));
    if (data.is_ho_geometry) {
      pipeline.push_back(new OpSetHOWeightsOnFace());
      pipeline.push_back(new OpGetHONormalsOnFace("MESH_NODE_POSITIONS"));
      pipeline.push_back(new OpCalculateHOCoords("MESH_NODE_POSITIONS"));
    }
    pipeline.push_back(new OpHOSetContravariantPiolaTransformOnFace3D(HDIV));
    // pipeline.push_back(new OpHOSetCovariantPiolaTransformOnFace3D(HDIV));
 
    pipeline.push_back(
        new OpCalculateVectorFieldValues<3>("U", commonDataPtr->mDispPtr));
    if (data.with_contact) {
      pipeline.push_back(
          new OpGetClosestRigidBodyData(mField, "U", commonDataPtr, false));
      pipeline.push_back(new OpCalculateHVecTensorTrace<3, BoundaryEleOp>(
          "SIGMA", commonDataPtr->contactTractionPtr));
    }
    MoFEMFunctionReturnHot(0);
  };
 
  auto add_boundary_ops_lhs = [&](auto &pipeline) {
    pipeline.push_back(new OpSetBc("U", true, boundaryMarker));
 
    // if (data.with_plasticity)
    //   pipeline.push_back(new OpSchurBeginBoundary("U", commonDataPtr));
 
    if (data.with_contact) {
      pipeline.push_back(
          new OpConstrainBoundaryLhs_dU("SIGMA", "U", commonDataPtr));
      pipeline.push_back(new OpConstrainBoundaryLhs_dTraction("SIGMA", "SIGMA",
                                                              commonDataPtr));
      // FIXME: precondition matrix
      pipeline.push_back(new OpSchurPreconditionMassContactBound(
          "SIGMA", "SIGMA", [](double, double, double) { return 1.; },
          commonDataPtr));
    }
 
    // pipeline.push_back(new OpSpringLhs(
    //     "U", "U",
 
    //     [this](double, double, double) { return (*cache).spring_stiffness; }
 
    //     ));
 
    if (data.is_rotating) {
      auto vol_side_fe_ptr = boost::make_shared<MoFEM::DomainSideEle>(mField);
 
      if (data.is_ho_geometry)
        CHKERR add_ho_ops(vol_side_fe_ptr->getOpPtrVector());
 
      vol_side_fe_ptr->getOpPtrVector().push_back(
          new OpVolumeSideAssembleRowData("U", commonDataPtr));
 
      pipeline.push_back(new OpRotatingFrameBoundaryLhs("U", "U", commonDataPtr,
                                                        vol_side_fe_ptr));
    }
    if (data.with_plasticity)
      pipeline.push_back(new OpSchurEndBoundary("U", commonDataPtr));
 
    pipeline.push_back(new OpUnSetBc("U"));
  };
 
  auto add_boundary_ops_rhs = [&](auto &pipeline) {
    pipeline.push_back(new OpSetBc("U", true, boundaryMarker));
    if (data.debug_info && data.with_contact)
      pipeline.push_back(new OpGetGaussPtsContactState("SIGMA", commonDataPtr));
 
    if (data.with_contact)
      pipeline.push_back(new OpConstrainBoundaryRhs("SIGMA", commonDataPtr));
    // pipeline.push_back(new OpSpringRhs(
    //     "U", commonDataPtr->mDispPtr,
    //     [this](double, double, double) { return (*cache).spring_stiffness;
    //     }));
 
    if (data.is_rotating) {
      auto my_vol_side_fe_ptr =
          boost::make_shared<MoFEM::DomainSideEle>(mField);
      my_vol_side_fe_ptr->getOpPtrVector().push_back(
          new OpCalculateVectorFieldGradient<3, 3>("U",
                                                   commonDataPtr->mGradPtr));
      pipeline.push_back(new OpRotatingFrameBoundaryRhs("U", commonDataPtr,
                                                        my_vol_side_fe_ptr));
    }
    pipeline.push_back(new OpUnSetBc("U"));
  };
 
  auto &tD = commonDataPtr->mtD;
  auto &tD_axi = commonDataPtr->mtD_Axiator;
  auto &tD_dev = commonDataPtr->mtD_Deviator;
 
  feLambdaLhs = boost::make_shared<DomainEle>(mField);
  feLambdaRhs = boost::make_shared<DomainEle>(mField);
  auto &lam_pipeline_lhs = feLambdaLhs->getOpPtrVector();
  auto &lam_pipeline_rhs = feLambdaRhs->getOpPtrVector();
 
  // boundary
  CHKERR add_boundary_base_ops(pipeline_mng->getOpBoundaryLhsPipeline());
  CHKERR add_boundary_base_ops(pipeline_mng->getOpBoundaryRhsPipeline());
  CHKERR add_disp_bc_ops();
 
  // pipeline for mechanical
  if (data.is_lambda_split) {
 
    CHKERR add_domain_base_ops(pipeline_mng->getOpDomainLhsPipeline());
    CHKERR add_domain_base_ops(pipeline_mng->getOpDomainRhsPipeline());
    CHKERR add_domain_stress_ops(pipeline_mng->getOpDomainLhsPipeline(),
                                 tD_dev);
    CHKERR add_domain_stress_ops(pipeline_mng->getOpDomainRhsPipeline(),
                                 tD_dev);
    CHKERR add_domain_ops_lhs(pipeline_mng->getOpDomainLhsPipeline(), tD_dev);
    CHKERR add_domain_ops_rhs(pipeline_mng->getOpDomainRhsPipeline(), tD_dev);
    CHKERR add_domain_ops_lhs_constraints(
        pipeline_mng->getOpDomainLhsPipeline(), tD_dev);
    CHKERR add_domain_ops_rhs_constraints(
        pipeline_mng->getOpDomainRhsPipeline());
 
    CHKERR add_domain_base_ops(lam_pipeline_lhs); //
    CHKERR add_domain_base_ops(lam_pipeline_rhs); //
    CHKERR add_domain_stress_ops(lam_pipeline_lhs, tD_axi);
    CHKERR add_domain_stress_ops(lam_pipeline_rhs, tD_axi);
    CHKERR add_domain_ops_lhs(lam_pipeline_lhs, tD_axi);
    CHKERR add_domain_ops_rhs(lam_pipeline_rhs, tD_axi);
 
  } else {
 
    CHKERR add_domain_base_ops(pipeline_mng->getOpDomainLhsPipeline());
    CHKERR add_domain_base_ops(pipeline_mng->getOpDomainRhsPipeline());
    CHKERR add_domain_stress_ops(pipeline_mng->getOpDomainLhsPipeline(), tD);
    CHKERR add_domain_stress_ops(pipeline_mng->getOpDomainRhsPipeline(), tD);
    CHKERR add_domain_ops_lhs(pipeline_mng->getOpDomainLhsPipeline(), tD);
    CHKERR add_domain_ops_rhs(pipeline_mng->getOpDomainRhsPipeline(), tD);
 
    // pipeline for constraints (only mu-part)
    CHKERR add_domain_stress_ops(pipeline_mng->getOpDomainLhsPipeline(),
                                 tD_dev);
    CHKERR add_domain_stress_ops(pipeline_mng->getOpDomainRhsPipeline(),
                                 tD_dev);
    CHKERR add_domain_ops_lhs_constraints(
        pipeline_mng->getOpDomainLhsPipeline(), tD_dev);
    CHKERR add_domain_ops_rhs_constraints(
        pipeline_mng->getOpDomainRhsPipeline());
  }
 
  // contact integration volume
  add_domain_ops_contact_lhs(pipeline_mng->getOpDomainLhsPipeline());
  add_domain_ops_contact_rhs(pipeline_mng->getOpDomainRhsPipeline());
 
  if (data.is_rotating && data.with_plasticity && false)
    add_skeleton_base_ops(pipeline_mng->getOpSkeletonRhsPipeline());
 
  // contact integration on boundary
  add_boundary_ops_lhs(pipeline_mng->getOpBoundaryLhsPipeline());
  add_boundary_ops_rhs(pipeline_mng->getOpBoundaryRhsPipeline());
 
  // setting integration rules
  auto integration_rule_bound = [](int, int, int approx_order) {
    return 2 * approx_order + 1;
    return 3 * approx_order;
  };
 
  auto integration_rule = [](int, int, int approx_order) {
    return 2 * (approx_order - 1);
    return 3 * approx_order;
  };
 
  auto integration_rule_lambda = [](int, int, int approx_order) { return 1; };
 
  CHKERR pipeline_mng->setDomainRhsIntegrationRule(integration_rule);
  CHKERR pipeline_mng->setDomainLhsIntegrationRule(integration_rule);
  CHKERR pipeline_mng->setBoundaryRhsIntegrationRule(integration_rule_bound);
  CHKERR pipeline_mng->setBoundaryLhsIntegrationRule(integration_rule_bound);
 
  feLambdaLhs->getRuleHook = integration_rule_lambda;
  feLambdaRhs->getRuleHook = integration_rule_lambda;
 
  if (data.with_contact) {
    // store data about the roller
    int def_val = 0;
    CHKERR mField.get_moab().tag_get_handle(
        "_ROLLER_NB", 1, MB_TYPE_INTEGER, commonDataPtr->rollerTag,
        MB_TAG_CREAT | MB_TAG_SPARSE, &def_val);
 
    data.update_roller = boost::make_shared<BoundaryEle>(mField);
    auto &pipeline = data.update_roller->getOpPtrVector();
    if (data.is_ho_geometry) {
      pipeline.push_back(new OpGetHONormalsOnFace("MESH_NODE_POSITIONS"));
      pipeline.push_back(new OpCalculateHOCoords("MESH_NODE_POSITIONS"));
    }
    pipeline.push_back(new OpHOSetContravariantPiolaTransformOnFace3D(HDIV));
 
    pipeline.push_back(
        new OpCalculateVectorFieldValues<3>("U", commonDataPtr->mDispPtr));
    pipeline.push_back(new OpCalculateHVecTensorTrace<3, BoundaryEleOp>(
        "SIGMA", commonDataPtr->contactTractionPtr));
    pipeline.push_back(
        new OpGetClosestRigidBodyData(mField, "U", commonDataPtr, true));
    if (data.debug_info)
      pipeline.push_back(
          new OpPostProcVertex("SIGMA", commonDataPtr, data.moab_debug));
    pipeline.push_back(new OpGetContactArea("SIGMA", commonDataPtr));
    data.update_roller->getRuleHook = integration_rule_bound;
  }
 
  auto create_reactions_element = [&]() {
    MoFEMFunctionBeginHot;
    Range &reaction_ents = commonDataPtr->reactionEnts;
 
    auto get_reactions_meshset_ents = [&]() {
      MoFEMFunctionBeginHot;
      for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(mField, BLOCKSET, it))
        if (it->getMeshsetId() == data.reaction_id) {
          Range ents;
          CHKERR mField.get_moab().get_entities_by_handle(it->meshset, ents,
                                                          true);
          CHKERR mField.get_moab().get_connectivity(ents, reaction_ents, true);
          CHKERR mField.get_moab().get_adjacencies(
              ents, 1, false, reaction_ents, moab::Interface::UNION);
          reaction_ents.merge(ents);
          break;
        }
 
      CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
          reaction_ents);
      int lents_size, ents_size;
      lents_size = reaction_ents.size();
      auto dm = simple->getDM();
      CHKERR MPIU_Allreduce(&lents_size, &ents_size, 1, MPIU_INT, MPIU_SUM,
                            PetscObjectComm((PetscObject)dm));
      if (!ents_size && data.reaction_id != -1)
        MOFEM_LOG_C("WORLD", Sev::warning,
                    "Warning: Provided meshset (id: %d) for calculating "
                    "reactions is EMPTY!",
                    data.reaction_id);
      // SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
      //          "Provide correct meshset id for calculating reactions "
      //          "(-reaction_id). ID = %d does not exist on input mesh.",
      //          data.reaction_id);
 
      MoFEMFunctionReturnHot(0);
    };
 
    CHKERR get_reactions_meshset_ents();
 
    double defs[] = {0, 0, 0};
    CHKERR mField.get_moab().tag_get_handle("_REACTION_TAG", 3, MB_TYPE_DOUBLE,
                                            commonDataPtr->reactionTag,
                                            MB_TAG_CREAT | MB_TAG_SPARSE, defs);
 
    // create calculate reactions element
    data.calc_reactions = boost::make_shared<DomainEle>(mField);
 
    std::vector<int> ghosts{0, 1, 2};
    commonDataPtr->reactionsVec = createSmartGhostVector(
        mField.get_comm(), (mField.get_comm_rank() ? 0 : 3), 3,
        (mField.get_comm_rank() ? 3 : 0), &*ghosts.begin());
 
    data.calc_reactions->getRuleHook = integration_rule;
 
    auto &pipeline = data.calc_reactions->getOpPtrVector();
    pipeline.push_back(
        new OpSetMaterialBlock("U", commonDataPtr, data.mapBlockTets));
    if (data.is_ho_geometry)
      CHKERR add_ho_ops(pipeline);
 
    pipeline.push_back(
        new OpCalculateVectorFieldGradient<3, 3>("U", commonDataPtr->mGradPtr));
    if (data.is_large_strain) {
      pipeline.push_back(new OpLogStrain("U", commonDataPtr));
    } else
      pipeline.push_back(new OpStrain("U", commonDataPtr));
 
    if (data.with_plasticity) {
      pipeline.push_back(new OpCalculateTensor2SymmetricFieldValues<3>(
          "EP", commonDataPtr->plasticStrainPtr));
      pipeline.push_back(
          new OpPlasticStress("U", commonDataPtr, commonDataPtr->mtD));
    } else
      pipeline.push_back(new OpStress("U", commonDataPtr, commonDataPtr->mtD));
 
    // FIXME:  include body forces
    // pipeline.push_back(new OpBodyForce("U", commonDataPtr, gravity));
    // if (data.is_rotating) {
    //   pipeline.push_back(new OpRotatingFrameRhs("U", commonDataPtr));
    //   pipeline.push_back(new OpCentrifugalForce2("U", commonDataPtr));
    // }
 
    if (data.is_large_strain) {
      if (data.is_neohooke)
        pipeline.push_back(
            new OpCalculateNeoHookeStressTangent("U", commonDataPtr));
      else
        pipeline.push_back(new OpCalculateLogStressTangent<false>(
            "U", commonDataPtr, commonDataPtr->mtD));
      pipeline.push_back(
          new OpInternalForceRhs<true, true>("U", commonDataPtr));
    } else
      pipeline.push_back(
          new OpInternalForceRhs<false, true>("U", commonDataPtr));
 
    MoFEMFunctionReturnHot(0);
  };
 
  if (data.calculate_reactions)
    CHKERR create_reactions_element();
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode ContactPlasticity::tsSolve() {
  MoFEMFunctionBegin;
 
  Simple *simple = mField.getInterface<Simple>();
  PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
  ISManager *is_manager = mField.getInterface<ISManager>();
 
  auto preProc = boost::make_shared<FePrePostProcess>(
      mField, (*cache).rollerDataVec, bcData, commonDataPtr);
 
  auto create_custom_ts = [&]() {
    auto set_dm_section = [&](auto dm) {
      MoFEMFunctionBegin;
      PetscSection section;
      CHKERR mField.getInterface<ISManager>()->sectionCreate(
          simple->getProblemName(), &section);
      CHKERR DMSetDefaultSection(dm, section);
      CHKERR DMSetDefaultGlobalSection(dm, section);
      CHKERR PetscSectionDestroy(&section);
      MoFEMFunctionReturn(0);
    };
    auto dm = simple->getDM();
 
    CHKERR set_dm_section(dm);
    boost::shared_ptr<FEMethod> null;
    preProc->methodsOpForMove = boost::shared_ptr<MethodForForceScaling>(
        new TimeForceScale(data.move_history, false));
 
    preProc->methodsOpForRollersDisp = boost::shared_ptr<MethodForForceScaling>(
        new TimeForceScale(data.contact_history, false));
    // here we use accelerogram as position data
    for (auto &roller : (*cache).rollerDataVec) {
      if (!roller.positionDataParamName.empty())
        preProc->methodsOpForRollersPosition.push_back(
            boost::shared_ptr<MethodForForceScaling>(
                new TimeAccelerogram(roller.positionDataParamName)));
    }
 
    // array<Range, 3> bc_ents;
    // for (auto &bdata : bcData)
    //   bc_ents[bdata.comp].merge(bdata.ents);
 
    // auto eraseRows = boost::make_shared<EraseRows>(
    //     mField, simple->getProblemName(), bc_ents[0], bc_ents[1],
    //     bc_ents[2]);
 
    if (true) {
 
      // Add element to calculate lhs of stiff part
      CHKERR DMMoFEMTSSetIJacobian(dm, simple->getDomainFEName(), null,
                                   dirichlet_bc_ptr, null);
      CHKERR DMMoFEMTSSetIJacobian(dm, simple->getDomainFEName(), null, preProc,
                                   null);
 
      auto postproc_bound_el = [&]() {
        MoFEMFunctionBeginHot;
        if (auto bmc_ptr = block_mat_container_ptr.lock()) {
          auto &bmc = *bmc_ptr;
          bmc.clear();
          // print schur complement after assembling the boundary contributions
          // if (commonDataPtr->STauTau)
          //     CHKERR MatView(commonDataPtr->STauTau,
          //     PETSC_VIEWER_DRAW_WORLD);
        }
        MoFEMFunctionReturnHot(0);
      };
 
      if (data.with_plasticity)
        pipeline_mng->getBoundaryLhsFE()->postProcessHook = postproc_bound_el;
 
      if (pipeline_mng->getBoundaryLhsFE())
        CHKERR DMMoFEMTSSetIJacobian(dm, simple->getBoundaryFEName(),
                                     pipeline_mng->getBoundaryLhsFE(), null,
                                     pipeline_mng->getBoundaryLhsFE());
 
      CHKERR DMMoFEMTSSetIJacobian(dm, simple->getDomainFEName(),
                                   pipeline_mng->getDomainLhsFE(), null, null);
      if (data.is_lambda_split)
        CHKERR DMMoFEMTSSetIJacobian(dm, simple->getDomainFEName(), feLambdaLhs,
                                     null, null);
 
      if (pipeline_mng->getSkeletonLhsFE())
        CHKERR DMMoFEMTSSetIJacobian(dm, simple->getSkeletonFEName(),
                                     pipeline_mng->getSkeletonLhsFE(), null,
                                     null);
 
      // CHKERR DMMoFEMTSSetIJacobian(dm, simple->getDomainFEName(), null, null,
      //                              dirichlet_bc_ptr);
      CHKERR DMMoFEMTSSetIJacobian(dm, simple->getDomainFEName(), null, null,
                                   preProc);
 
      if (pipeline_mng->getDomainRhsFE()) {
 
        // add dirichlet boundary conditions
        CHKERR DMMoFEMTSSetIFunction(dm, simple->getDomainFEName(), null,
                                     preProc, null);
        if (data.is_lambda_split)
          CHKERR DMMoFEMTSSetIFunction(dm, simple->getDomainFEName(),
                                       feLambdaRhs, null, null);
        CHKERR DMMoFEMTSSetIFunction(dm, simple->getDomainFEName(), null,
                                     dirichlet_bc_ptr, null);
 
        CHKERR DMMoFEMTSSetIFunction(dm, simple->getDomainFEName(),
                                     pipeline_mng->getDomainRhsFE(), null,
                                     null);
 
        if (pipeline_mng->getBoundaryRhsFE())
          CHKERR DMMoFEMTSSetIFunction(dm, simple->getBoundaryFEName(),
                                       pipeline_mng->getBoundaryRhsFE(), null,
                                       null);
        if (pipeline_mng->getSkeletonRhsFE())
          CHKERR DMMoFEMTSSetIFunction(dm, simple->getSkeletonFEName(),
                                       pipeline_mng->getSkeletonRhsFE(), null,
                                       null);
        // Add surface forces
        for (auto fit = neumann_forces.begin(); fit != neumann_forces.end();
             fit++) {
          CHKERR DMMoFEMTSSetIFunction(dm, fit->first.c_str(),
                                       &fit->second->getLoopFe(), NULL, NULL);
        }
 
        // Add edge forces
        for (auto fit = edge_forces.begin(); fit != edge_forces.end(); fit++) {
          CHKERR DMMoFEMTSSetIFunction(dm, fit->first.c_str(),
                                       &fit->second->getLoopFe(), NULL, NULL);
        }
 
        // Add nodal forces
        for (auto fit = nodal_forces.begin(); fit != nodal_forces.end();
             fit++) {
          CHKERR DMMoFEMTSSetIFunction(dm, fit->first.c_str(),
                                       &fit->second->getLoopFe(), NULL, NULL);
        }
 
        // CHKERR DMMoFEMTSSetIFunction(dm, simple->getDomainFEName(), null,
        // null,
        //                              dirichlet_bc_ptr);
        CHKERR DMMoFEMTSSetIFunction(dm, simple->getDomainFEName(), null, null,
                                     preProc);
      }
    }
 
    auto print_active_set = [&](std::array<int, 2> &lgauss_pts, string name,
                                MoFEM::Interface &mField) {
      MoFEMFunctionBeginHot;
      std::vector<int> gauss_pts(2, 0);
      auto dm = mField.getInterface<Simple>()->getDM();
      CHKERR MPIU_Allreduce(lgauss_pts.data(), gauss_pts.data(), 2, MPIU_INT,
                            MPIU_SUM, PetscObjectComm((PetscObject)dm));
 
      MOFEM_LOG_C("WORLD", Sev::inform, "\t \t Active %s gauss pts: %d / %d",
                  name.c_str(), (int)gauss_pts[0], (int)gauss_pts[1]);
      lgauss_pts[0] = lgauss_pts[1] = 0;
      MoFEMFunctionReturnHot(0);
    };
 
    auto postproc_dom = [&]() {
      MoFEMFunctionBeginHot;
      CHKERR print_active_set(commonDataPtr->stateArrayPlast, "plastic",
                              mField);
      MoFEMFunctionReturnHot(0);
    };
    auto postproc_bound = [&]() {
      MoFEMFunctionBeginHot;
      CHKERR print_active_set(commonDataPtr->stateArrayCont, "contact", mField);
      MoFEMFunctionReturnHot(0);
    };
 
    if (data.debug_info) {
      if (data.with_plasticity)
        pipeline_mng->getDomainRhsFE()->postProcessHook = postproc_dom;
 
      if (data.with_contact)
        pipeline_mng->getBoundaryRhsFE()->postProcessHook = postproc_bound;
    }
 
    auto ts = MoFEM::createTS(mField.get_comm());
    CHKERR TSSetDM(ts, dm);
    return ts;
  };
 
  auto solver = create_custom_ts();
 
  CHKERR TSSetExactFinalTime(solver, TS_EXACTFINALTIME_MATCHSTEP);
 
  auto dm = simple->getDM();
  auto D = smartCreateDMVector(dm);
 
  CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR TSSetFromOptions(solver);
  CHKERR TSSetSolution(solver, D);
 
  CHKERR TSSetUp(solver);
 
  // constexpr bool run_fieldsplit_for_boundary = true;
  if (true) {
    SNES snes;
    CHKERR TSGetSNES(solver, &snes);
    KSP ksp;
    CHKERR SNESGetKSP(snes, &ksp);
    PC pC;
    CHKERR KSPGetPC(ksp, &pC);
    PetscBool is_pcfs = PETSC_FALSE;
    PetscObjectTypeCompare((PetscObject)pC, PCFIELDSPLIT, &is_pcfs);
    SmartPetscObj<AO> sigma_ao; // contact ordering
 
    auto make_is_field_map = [&]() {
      PetscSection section;
      CHKERR DMGetDefaultSection(dm, &section);
 
      int num_fields;
      CHKERR PetscSectionGetNumFields(section, &num_fields);
 
      const MoFEM::Problem *prb_ptr;
      CHKERR DMMoFEMGetProblemPtr(dm, &prb_ptr);
 
      map<std::string, SmartPetscObj<IS>> is_map;
      for (int ff = 0; ff != num_fields; ff++) {
 
        const char *field_name;
        CHKERR PetscSectionGetFieldName(section, ff, &field_name);
 
        CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
            prb_ptr->getName(), ROW, field_name, 0, MAX_DOFS_ON_ENTITY,
            is_map[field_name]);
      };
      return is_map;
    };
    auto is_map = make_is_field_map();
 
    auto set_fieldsplit_on_bc = [&](PC &bc_pc, string name_prb) {
      MoFEMFunctionBeginHot;
 
      PetscBool is_bc_field_split;
      PetscObjectTypeCompare((PetscObject)bc_pc, PCFIELDSPLIT,
                             &is_bc_field_split);
 
      // check if we have bc ents on any proc
      int l_has_bcs_ents = allBoundaryEnts->size();
      int g_has_bcs_ents;
      CHKERR MPIU_Allreduce(&l_has_bcs_ents, &g_has_bcs_ents, 1, MPIU_INT,
                            MPIU_SUM, PetscObjectComm((PetscObject)dm));
 
      if (is_bc_field_split && g_has_bcs_ents) {
        MOFEM_LOG_C("WORLD", Sev::inform,
                    "Applying fieldsplit preconditioner for %d "
                    "boundary entities.",
                    g_has_bcs_ents);
 
        Range nodes;
        CHKERR mField.get_moab().get_connectivity(*allBoundaryEnts, nodes,
                                                  true);
        allBoundaryEnts->merge(nodes);
        CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
            *allBoundaryEnts);
 
        SmartPetscObj<IS> is_bc;
        CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
            name_prb, ROW, "U", 0, 3, is_bc, allBoundaryEnts.get());
 
        CHKERR ISSort(is_bc);
        CHKERR PCFieldSplitSetIS(bc_pc, PETSC_NULL,
                                 is_bc); // boundary block
        CHKERR PCSetUp(bc_pc);
      } else {
        CHKERR PCSetFromOptions(bc_pc);
        CHKERR PCSetType(bc_pc, PCLU);
        CHKERR PCSetUp(bc_pc);
      }
 
      MoFEMFunctionReturnHot(0);
    };
 
    auto set_global_mat_and_vec = [&]() {
      MoFEMFunctionBeginHot;
      auto fe = mField.getInterface<PipelineManager>()->getDomainLhsFE();
      CHKERR DMCreateMatrix_MoFEM(dm, &fe->ts_B);
      CHKERR TSSetIFunction(solver, fe->ts_F, PETSC_NULL, PETSC_NULL);
      CHKERR TSSetIJacobian(solver, fe->ts_B, fe->ts_B, PETSC_NULL, PETSC_NULL);
      CHKERR KSPSetOperators(ksp, fe->ts_B, fe->ts_B);
      MoFEMFunctionReturnHot(0);
    };
 
    // Setup fieldsplit (block) solver - optional: yes/no
    if (is_pcfs == PETSC_TRUE && !data.is_fieldsplit_bc_only) {
 
      CHKERR set_global_mat_and_vec();
 
      PetscBool is_l0_field_split;              // level0 fieldsplit
      PetscBool is_l1_field_split = PETSC_TRUE; // level1 fieldsplit
 
      PetscObjectTypeCompare((PetscObject)pC, PCFIELDSPLIT, &is_l0_field_split);
    
      auto set_fieldsplit_contact = [&]() {
        MoFEMFunctionBeginHot;
 
        const MoFEM::Problem *fs_sig_ptr;
        CHKERR DMMoFEMGetProblemPtr(dmContact, &fs_sig_ptr);
        if (auto sig_data = fs_sig_ptr->subProblemData) {
 
          is_map["E_IS_SIG"] = sig_data->rowIs;
          sigma_ao = sig_data->getSmartRowMap();
 
          CHKERR PCFieldSplitSetIS(pC, NULL, is_map["SIGMA"]);
          CHKERR PCFieldSplitSetIS(pC, NULL, is_map["E_IS_SIG"]);
          //  CHKERR PCFieldSplitSetType(pC, PC_COMPOSITE_ADDITIVE);
          //  CHKERR PCFieldSplitSetType(pC,
          //  PC_COMPOSITE_MULTIPLICATIVE);
          CHKERR PCSetUp(pC);
 
          SmartPetscObj<Mat> mat_test;
          CHKERR DMCreateMatrix_MoFEM(dmContact, mat_test);
 
          PetscInt n;
          KSP *ct_ksp;
          CHKERR PCFieldSplitGetSubKSP(pC, &n, &ct_ksp);
          PC l1_pc;
          CHKERR KSPGetPC(ct_ksp[1], &l1_pc);
          PetscObjectTypeCompare((PetscObject)l1_pc, PCFIELDSPLIT,
                                 &is_l1_field_split);
          // important!
          // in case of contact analysis we swap the first preconditioner
          pC = l1_pc;
          CHKERR PetscFree(ct_ksp);
        }
        MoFEMFunctionReturnHot(0);
      };
 
      if (data.with_contact && is_l0_field_split)
        CHKERR set_fieldsplit_contact();
 
      PetscObjectTypeCompare((PetscObject)pC, PCFIELDSPLIT, &is_l1_field_split);
      if (data.with_plasticity && is_l0_field_split && is_l1_field_split) {
 
        const MoFEM::Problem *schur_ep_ptr;
        CHKERR DMMoFEMGetProblemPtr(dmEP, &schur_ep_ptr);
        if (auto ep_data = schur_ep_ptr->subProblemData) {
 
          is_map["E_IS_EP"] = ep_data->rowIs;
          // CHKERR ISView(is_map["EP"], PETSC_VIEWER_STDOUT_SELF);
          if (sigma_ao) 
            CHKERR AOApplicationToPetscIS(sigma_ao, is_map["EP"]);
          // CHKERR ISSort(is_map["EP"]);
 
          CHKERR PCFieldSplitSetIS(pC, NULL, is_map["EP"]);
          CHKERR PCFieldSplitSetIS(pC, NULL, is_map["E_IS_EP"]);
 
          PCCompositeType pc_type;
          CHKERR PCFieldSplitGetType(pC, &pc_type);
 
          if (pc_type == PC_COMPOSITE_SCHUR) {
            CHKERR DMCreateMatrix_MoFEM(dmEP, commonDataPtr->SEpEp);
            commonDataPtr->aoSEpEp = ep_data->getSmartRowMap();
            if (sigma_ao)
              commonDataPtr->aoSEpEp = sigma_ao;
 
            commonDataPtr->blockMatContainerPtr =
                boost::make_shared<BlockMatContainer>();
            MoFEM::block_mat_container_ptr =
                commonDataPtr->blockMatContainerPtr;
 
            CHKERR PCFieldSplitSetSchurPre(pC, PC_FIELDSPLIT_SCHUR_PRE_USER,
                                           commonDataPtr->SEpEp);
 
            CHKERR PCSetUp(pC);
          }
 
          PetscInt n;
          KSP *tt_ksp;
          CHKERR PCFieldSplitGetSubKSP(pC, &n, &tt_ksp);
          PC tau_pc;
          CHKERR KSPGetPC(tt_ksp[1], &tau_pc);
          PetscBool is_tau_field_split;
          PetscObjectTypeCompare((PetscObject)tau_pc, PCFIELDSPLIT,
                                 &is_tau_field_split);
 
          if (is_tau_field_split) {
            const MoFEM::Problem *schur_tau_ptr;
            CHKERR DMMoFEMGetProblemPtr(dmTAU, &schur_tau_ptr);
            if (auto tau_data = schur_tau_ptr->subProblemData) {
 
              // CHKERR tau_data->getRowIs(is_map["E_IS_TAU"]);
              is_map["E_IS_TAU"] = tau_data->rowIs;
 
              AO ep_ao;
              CHKERR ep_data->getRowMap(&ep_ao);
 
              if (sigma_ao)
                CHKERR AOApplicationToPetscIS(sigma_ao, is_map["TAU"]);
              CHKERR AOApplicationToPetscIS(ep_ao, is_map["TAU"]);
 
              CHKERR PCFieldSplitSetIS(tau_pc, NULL, is_map["TAU"]);
              CHKERR PCFieldSplitSetIS(tau_pc, NULL, is_map["E_IS_TAU"]);
 
              CHKERR PCFieldSplitGetType(tau_pc, &pc_type);
              if (pc_type == PC_COMPOSITE_SCHUR) {
                CHKERR DMCreateMatrix_MoFEM(dmTAU, commonDataPtr->STauTau);
                commonDataPtr->aoSTauTau = tau_data->getSmartRowMap();
 
                CHKERR PCFieldSplitSetSchurPre(tau_pc,
                                               PC_FIELDSPLIT_SCHUR_PRE_USER,
                                               commonDataPtr->STauTau);
              }
              CHKERR PCSetUp(tau_pc);
 
              PetscInt bc_n;
              KSP *bc_ksp;
              CHKERR PCFieldSplitGetSubKSP(tau_pc, &bc_n, &bc_ksp);
              PC bc_pc;
              CHKERR KSPGetPC(bc_ksp[1], &bc_pc);
 
              CHKERR set_fieldsplit_on_bc(bc_pc, schur_tau_ptr->getName());
              CHKERR PetscFree(bc_ksp);
            }
          }
          CHKERR PetscFree(tt_ksp);
        }
      }
    } else if (is_pcfs && data.is_fieldsplit_bc_only) {
 
      char mumps_options[] = "-fieldsplit_0_mat_mumps_icntl_14 1200 "
                             "-fieldsplit_0_mat_mumps_icntl_24 1 "
                             "-fieldsplit_0_mat_mumps_icntl_13 1 "
                             "-fieldsplit_1_mat_mumps_icntl_14 1200 "
                             "-fieldsplit_1_mat_mumps_icntl_24 1 "
                             "-fieldsplit_1_mat_mumps_icntl_13 1";
      CHKERR PetscOptionsInsertString(NULL, mumps_options);
      CHKERR set_global_mat_and_vec();
      CHKERR set_fieldsplit_on_bc(pC, simple->getProblemName());
 
    } else {
      SNES snes;
      CHKERR TSGetSNES(solver, &snes);
      KSP ksp;
      CHKERR SNESGetKSP(snes, &ksp);
      PC pCc;
      CHKERR KSPGetPC(ksp, &pCc);
      CHKERR PCSetFromOptions(pCc);
      CHKERR PCSetType(pCc, PCLU);
    }
  }
 
  auto set_section_monitor = [&]() {
    MoFEMFunctionBegin;
    SNES snes;
    CHKERR TSGetSNES(solver, &snes);
    PetscViewerAndFormat *vf;
    CHKERR PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,
                                      PETSC_VIEWER_DEFAULT, &vf);
    CHKERR SNESMonitorSet(
        snes,
        (MoFEMErrorCode(*)(SNES, PetscInt, PetscReal, void *))SNESMonitorFields,
        vf, (MoFEMErrorCode(*)(void **))PetscViewerAndFormatDestroy);
    MoFEMFunctionReturn(0);
  };
 
  auto create_post_process_element = [&]() {
    MoFEMFunctionBegin;
    postProcFe = boost::make_shared<PostProcVolumeOnRefinedMesh>(mField);
    postProcFe->generateReferenceElementMesh();
 
    postProcFe->getOpPtrVector().push_back(
        new OpSetMaterialBlock("U", commonDataPtr, data.mapBlockTets));
 
    postProcFe->getOpPtrVector().push_back(
        new OpCalculateVectorFieldGradient<3, 3>("U", commonDataPtr->mGradPtr));
    if (data.is_large_strain) {
      postProcFe->getOpPtrVector().push_back(
          new OpLogStrain("U", commonDataPtr));
    } else
      postProcFe->getOpPtrVector().push_back(new OpStrain("U", commonDataPtr));
    if (data.with_plasticity) {
      postProcFe->getOpPtrVector().push_back(new OpCalculateScalarFieldValues(
          "TAU", commonDataPtr->plasticTauPtr));
      postProcFe->getOpPtrVector().push_back(
          new OpCalculateTensor2SymmetricFieldValues<3>(
              "EP", commonDataPtr->plasticStrainPtr));
      postProcFe->getOpPtrVector().push_back(
          new OpPlasticStress("U", commonDataPtr, commonDataPtr->mtD));
      postProcFe->getOpPtrVector().push_back(
          new OpCalculatePlasticSurfaceAndFlow("U", commonDataPtr));
    } else
      postProcFe->getOpPtrVector().push_back(
          new OpStress("U", commonDataPtr, commonDataPtr->mtD));
 
    if (data.is_large_strain) {
      if (data.is_neohooke)
        postProcFe->getOpPtrVector().push_back(
            new OpCalculateNeoHookeStressTangent("U", commonDataPtr));
      else
        postProcFe->getOpPtrVector().push_back(
            new OpCalculateLogStressTangent<false>("U", commonDataPtr,
                                                   commonDataPtr->mtD));
    }
    if (data.with_contact) {
      postProcFe->getOpPtrVector().push_back(
          new OpCalculateHVecTensorDivergence<3, 3>(
              "SIGMA", commonDataPtr->contactStressDivergencePtr));
      postProcFe->getOpPtrVector().push_back(
          new OpCalculateHVecTensorField<3, 3>(
              "SIGMA", commonDataPtr->contactStressPtr));
 
      postProcFe->getOpPtrVector().push_back(
          new OpPostProcContact("SIGMA", postProcFe->postProcMesh,
                                postProcFe->mapGaussPts, commonDataPtr));
    }
    if (data.is_large_strain)
      postProcFe->getOpPtrVector().push_back(
          new OpPostProcElastic<true>("U", postProcFe->postProcMesh,
                                      postProcFe->mapGaussPts, commonDataPtr));
    else
      postProcFe->getOpPtrVector().push_back(
          new OpPostProcElastic<false>("U", postProcFe->postProcMesh,
                                       postProcFe->mapGaussPts, commonDataPtr));
 
    if (data.with_plasticity)
      postProcFe->getOpPtrVector().push_back(
          new OpPostProcPlastic("U", postProcFe->postProcMesh,
                                postProcFe->mapGaussPts, commonDataPtr));
 
    if (data.calculate_reactions)
      // FIXME: use addFieldValuesPostProc("U", "REACTIONS", vec_rhs); instead
      postProcFe->getOpPtrVector().push_back(
          new OpSaveReactionForces(mField, "U", postProcFe->postProcMesh,
                                   postProcFe->mapGaussPts, commonDataPtr));
    postProcFe->addFieldValuesPostProc("U", "DISPLACEMENT");
    if (data.is_ho_geometry)
      postProcFe->addFieldValuesPostProc("MESH_NODE_POSITIONS");
    MoFEMFunctionReturn(0);
  };
 
  auto scatter_create = [&](auto coeff) {
    SmartPetscObj<IS> is;
    CHKERR is_manager->isCreateProblemFieldAndRank(simple->getProblemName(),
                                                   ROW, "U", coeff, coeff, is);
    int loc_size;
    CHKERR ISGetLocalSize(is, &loc_size);
    Vec v;
    CHKERR VecCreateMPI(mField.get_comm(), loc_size, PETSC_DETERMINE, &v);
    VecScatter scatter;
    CHKERR VecScatterCreate(D, is, v, PETSC_NULL, &scatter);
    return std::make_tuple(SmartPetscObj<Vec>(v),
                           SmartPetscObj<VecScatter>(scatter));
  };
 
  auto set_time_monitor = [&]() {
    MoFEMFunctionBegin;
    boost::shared_ptr<MMonitor> monitor_ptr(
        new MMonitor(dm, mField, commonDataPtr, postProcFe, uXScatter,
                     uYScatter, uZScatter));
    boost::shared_ptr<ForcesAndSourcesCore> null;
    CHKERR DMMoFEMTSSetMonitor(dm, solver, simple->getDomainFEName(),
                               monitor_ptr, null, null);
    MoFEMFunctionReturn(0);
  };
 
  CHKERR set_section_monitor();
  CHKERR create_post_process_element();
 
  if (data.is_restart) {
 
    CHKERR PetscPrintf(PETSC_COMM_WORLD,
                       "Reading vector in binary from %s file...\n",
                       data.restart_file_str.c_str());
    PetscViewer viewer;
    CHKERR PetscViewerBinaryOpen(PETSC_COMM_WORLD,
                                 data.restart_file_str.c_str(), FILE_MODE_READ,
                                 &viewer);
    CHKERR VecLoad(D, viewer);
 
    CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
 
    typedef boost::tokenizer<boost::char_separator<char>> Tokenizer;
    string s = data.restart_file_str;
    Tokenizer tok{s, boost::char_separator<char>("_")};
    auto it = tok.begin();
    const int restart_step = std::stoi((++it)->c_str());
    data.restart_step = restart_step;
    string restart_tt = *(++it);
    restart_tt.erase(restart_tt.length() - 4);
    const double restart_time = std::atof(restart_tt.c_str());
    CHKERR TSSetStepNumber(solver, restart_step);
    CHKERR TSSetTime(solver, restart_time);
  }
 
  uXScatter = scatter_create(0);
  uYScatter = scatter_create(1);
  uZScatter = scatter_create(2);
  CHKERR set_time_monitor();
 
  // Add iteration post-proc for debugging
  // boost::shared_ptr<FEMethod> null_fe;
  // SNES snes;
  // CHKERR TSGetSNES(solver, &snes);
  // auto write_fe = boost::make_shared<FEMethod>();
  // write_fe->postProcessHook = [&]() {
  //   MoFEMFunctionBegin;
  //   int iter;
  //   CHKERR SNESGetIterationNumber(snes, &iter);
  //   CHKERR postProcFe->writeFile(
  //       "it_out_plastic_" + boost::lexical_cast<std::string>(iter) + ".h5m");
  //   MoFEMFunctionReturn(0);
  // };
  // CHKERR DMMoFEMTSSetIJacobian(dm, simple->getDomainFEName(), postProcFe,
  //                              null_fe, write_fe);
 
  CHKERR TSSolve(solver, D);
 
  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode ContactPlasticity::postProcess() { return 0; }
 
MoFEMErrorCode ContactPlasticity::checkResults() { return 0; }
 
Range ContactPlasticity::getEntsOnMeshSkin() {
  Range body_ents;
  auto meshset = mField.getInterface<Simple>()->getMeshSet();
  CHKERR mField.get_moab().get_entities_by_dimension(meshset, 3, body_ents);
  // cout << "body_ents for skin: " << body_ents.size() << endl;
  Skinner skin(&mField.get_moab());
  Range skin_tris;
  CHKERR skin.find_skin(0, body_ents, false, skin_tris);
  Range boundary_ents;
  ParallelComm *pcomm =
      ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
  if (pcomm == NULL) {
    pcomm = new ParallelComm(&mField.get_moab(), mField.get_comm());
  }
  CHKERR pcomm->filter_pstatus(skin_tris, PSTATUS_SHARED | PSTATUS_MULTISHARED,
                               PSTATUS_NOT, -1, &boundary_ents);
 
  return boundary_ents;
};
 
static char help[] = "...\n\n";
 
int main(int argc, char *argv[]) {
 
#include <DefaultParams.hpp>
 
  MoFEM::Core::Initialize(&argc, &argv, param_file.c_str(), help);
 
  try {
 
    DMType dm_name = "DMMOFEM";
    CHKERR DMRegister_MoFEM(dm_name);
 
    moab::Core mb_instance;              ///< mesh database
    moab::Interface &moab = mb_instance; ///< mesh database interface
 
    MoFEM::Core core(moab);           ///< finite element database
    MoFEM::Interface &m_field = core; ///< finite element database insterface
 
    Simple *simple = m_field.getInterface<Simple>();
    simple->getProblemName() = "Multifield plasticity";
    CHKERR simple->getOptions();
 
    ContactPlasticity ex(m_field);
    CHKERR ContactPlasticity::getAnalysisParameters();
 
    PetscBool is_partitioned = PETSC_FALSE;
    CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-is_partitioned",
                               &is_partitioned, PETSC_NULL);
    if (is_partitioned) {
      CHKERR simple->loadFile();
    } else {
      if (m_field.get_comm_size() > 1 && true) {
        CHKERR simple->loadFile("", "");
        CHKERR ContactPlasticity::myMeshPartition(m_field);
      } else
        CHKERR simple->loadFile("");
    }
 
    // Add meshsets with material and boundary conditions
    CHKERR m_field.getInterface<MeshsetsManager>()->setMeshsetFromFile();
 
    CHKERR ex.loadMeshBlockSets();
    CHKERR cache->getOptionsFromCommandLine();
    if (ex.data.scale_params)
      for (auto &mit : mat_blocks)
        CHKERR mit.second.scaleParameters();
 
    CHKERR ex.runProblem();
  }
  CATCH_ERRORS;
 
  CHKERR MoFEM::Core::Finalize();
}