ContactPlasticity Struct Reference

Collaboration diagram for ContactPlasticity:

Classes
struct	MMonitor
struct	ProblemData

Public Member Functions
	ContactPlasticity (MoFEM::Interface &m_field)
MoFEMErrorCode	runProblem ()
MoFEMErrorCode	loadMeshBlockSets ()

Static Public Member Functions
static MoFEMErrorCode	getAnalysisParameters ()
static MoFEMErrorCode	myMeshPartition (MoFEM::Interface &m_field)

Public Attributes
FieldApproximationBase	base

Static Public Attributes
static ProblemData	data

Private Member Functions
MoFEMErrorCode	setUP ()
MoFEMErrorCode	createCommonData ()
MoFEMErrorCode	OPs ()
MoFEMErrorCode	bC ()
MoFEMErrorCode	tsSolve ()
MoFEMErrorCode	postProcess ()
MoFEMErrorCode	checkResults ()
Range	getEntsOnMeshSkin ()

Private Attributes
MoFEM::Interface &	mField
MatrixDouble	invJac
MatrixDouble	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
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::ptr_vector< DataFromMove >	bcData
boost::shared_ptr< DomainSideEle >	volSideFe
boost::shared_ptr< DomainEle >	feLambdaRhs
boost::shared_ptr< DomainEle >	feLambdaLhs
SmartPetscObj< DM >	dmEP
SmartPetscObj< DM >	dmTAU
SmartPetscObj< DM >	dmContact

Detailed Description

Definition at line 107 of file multifield_plasticity.cpp.

Constructor & Destructor Documentation

ContactPlasticity()

ContactPlasticity::ContactPlasticity ( MoFEM::Interface &  m_field )

inline

Definition at line 109 of file multifield_plasticity.cpp.

: mField(m_field) {}

Member Function Documentation

bC()

MoFEMErrorCode ContactPlasticity::bC ( )

private

Definition at line 1116 of file multifield_plasticity.cpp.

                                     {
  MoFEMFunctionBegin;
  auto bc_mng = mField.getInterface<BcManager>();
  auto prb_mng = mField.getInterface<ProblemsManager>();
  auto simple = mField.getInterface<Simple>();
 
  auto get_ents = [&](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;
  };
 
  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "FIX_X",
                                           "U", 0, 0);
  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "FIX_Y",
                                           "U", 1, 1);
  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "FIX_Z",
                                           "U", 2, 2);
  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "FIX_ALL",
                                           "U", 0, 2);
                                         
  if (data.with_contact) {
    Range boundary_ents;
    for (auto &bc : mField.getInterface<BcManager>()->getBcMapByBlockName())
      if (bc_mng->checkBlock(bc, "FIX_"))
        boundary_ents.merge(*bc.second->getBcEntsPtr());
    boundary_ents.merge(get_ents("NO_CONTACT"));
    CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
        boundary_ents);
    CHKERR mField.getInterface<ProblemsManager>()->removeDofsOnEntities(
        mField.getInterface<Simple>()->getProblemName(), "SIGMA", boundary_ents,
        0, 2);
  }
 
  std::vector<DataFromBc> bcData;
  dirichlet_bc_ptr = boost::make_shared<DirichletDisplacementRemoveDofsBc>(
      mField, "U", simple->getProblemName(), "DISPLACEMENT", true);
  dirichlet_bc_ptr->iNitialize();
 
  auto dm = simple->getDM();
  // FIXME: contact with fieldsplit does not work in parallel, requires further
  // investigation
  // TODO: will come back into this in the future
  if (data.with_contact && data.with_plasticity && false) {
    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);
  }
  // CHKERR DMoFEMPreProcessFiniteElements(dm, dirichlet_bc_ptr.get());
 
  // 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);
}

checkResults()

MoFEMErrorCode ContactPlasticity::checkResults ( )

private

Definition at line 2581 of file multifield_plasticity.cpp.

{ return 0; }

createCommonData()

MoFEMErrorCode ContactPlasticity::createCommonData ( )

private

Definition at line 1001 of file multifield_plasticity.cpp.

                                                   {
  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);
}

getAnalysisParameters()

MoFEMErrorCode ContactPlasticity::getAnalysisParameters ( )

static

Definition at line 602 of file multifield_plasticity.cpp.

                                                        {
  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_reduced_integration",
      "use axiator split approach, with reduced integration", "",
      data.is_reduced_integration, &data.is_reduced_integration, PETSC_NULL);
  CHKERR PetscOptionsBool("-test_user_base_tau", "test_user_base_for_tau", "",
                          data.test_user_base_tau, &data.test_user_base_tau,
                          PETSC_NULL);
  CHKERR PetscOptionsBool("-is_alm",
                          "use Augmented Lagrangian method for enforcing the "
                          "KKT condtions for plasticity",
                          "", data.is_alm, &data.is_alm, 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);
}

getEntsOnMeshSkin()

Range ContactPlasticity::getEntsOnMeshSkin ( )

private

Definition at line 2583 of file multifield_plasticity.cpp.

                                           {
  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;
};

loadMeshBlockSets()

MoFEMErrorCode ContactPlasticity::loadMeshBlockSets

(

)

Definition at line 793 of file multifield_plasticity.cpp.

                                                    {
  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);
}

myMeshPartition()

MoFEMErrorCode ContactPlasticity::myMeshPartition ( MoFEM::Interface &  m_field )

static

Definition at line 704 of file multifield_plasticity.cpp.

                                                                         {
  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);
}

OPs()

MoFEMErrorCode ContactPlasticity::OPs ( )

private

Definition at line 1217 of file multifield_plasticity.cpp.

                                      {
  MoFEMFunctionBegin;
  PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
  Simple *simple = mField.getInterface<Simple>();
  auto problem_manager = mField.getInterface<ProblemsManager>();
  auto bc_mng = mField.getInterface<BcManager>();
  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;
  };
 
  CHKERR bc_mng->pushMarkDOFsOnEntities(simple->getProblemName(), "MOVE_X", "U",
                                        0, 0);
  CHKERR bc_mng->pushMarkDOFsOnEntities(simple->getProblemName(), "MOVE_Y", "U",
                                        1, 1);
  CHKERR bc_mng->pushMarkDOFsOnEntities(simple->getProblemName(), "MOVE_Z", "U",
                                        2, 2);
 
  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();
  boundaryMarker =
      bc_mng->getMergedBlocksMarker(vector<string>{"MOVE_", "ROTATE_"});
 
  auto add_disp_bc_ops = [&]() {
    MoFEMFunctionBeginHot;
    std::map<char, size_t> xyz_map = {{'X', 0}, {'Y', 1}, {'Z', 2}};
 
    if (data.with_plasticity)
      pipeline_mng->getOpBoundaryLhsPipeline().push_back(
          new OpSchurBeginBoundary("U", commonDataPtr));
 
    for (auto &bc : mField.getInterface<BcManager>()->getBcMapByBlockName()) {
      if (bc_mng->checkBlock(bc, "MOVE_")) {
        
        int idx = xyz_map.at(bc.first[(bc.first.find("MOVE_") + 5)]);
 
        auto disp_vec = bc.second->bcAttributes;
        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(idx) = disp_vec[0];
 
        bcData.push_back(new DataFromMove(disp, *bc.second->getBcEntsPtr()));
        // rhs
        pipeline_mng->getOpBoundaryRhsPipeline().push_back(
            new OpSetBc("U", false, bc.second->getBcMarkersPtr()));
        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, bc.second->getBcMarkersPtr()));
        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"));
      }
    }
 
    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<3>(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);
 
  if (data.test_user_base_tau) {
    feLambdaLhs->getUserPolynomialBase() =
        boost::shared_ptr<BaseFunction>(new BubbleTauPolynomialBase());
    feLambdaRhs->getUserPolynomialBase() =
        boost::shared_ptr<BaseFunction>(new BubbleTauPolynomialBase());
  }
 
  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_reduced_integration) {
 
    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 = [this](int, int, int approx_order) {
    if (base == DEMKOWICZ_JACOBI_BASE)
      return 2 * approx_order - 1;
    // return 3 * approx_order;
    return 2 * (approx_order - 1);
  };
 
  auto integration_rule_lambda = [](int, int, int approx_order) {
    // return 3 * 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 (-reaction_id: %d) for calculating "
            "reactions is EMPTY! Calculating reactions for all FIX_ blocks.",
            data.reaction_id);
        for (auto &bc : mField.getInterface<BcManager>()->getBcMapByBlockName())
          if (bc_mng->checkBlock(bc, "FIX_"))
            reaction_ents.merge(*bc.second->getBcEntsPtr());
        CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
            reaction_ents);
      }
 
      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);
}

postProcess()

MoFEMErrorCode ContactPlasticity::postProcess ( )

private

Definition at line 2579 of file multifield_plasticity.cpp.

{ return 0; }

runProblem()

MoFEMErrorCode ContactPlasticity::runProblem

(

)

Definition at line 591 of file multifield_plasticity.cpp.

                                             {
  MoFEMFunctionBegin;
  CHKERR setUP();
  CHKERR createCommonData();
  CHKERR bC();
  CHKERR OPs();
  CHKERR tsSolve();
  CHKERR postProcess();
  CHKERR checkResults();
  MoFEMFunctionReturn(0);
}

setUP()

MoFEMErrorCode ContactPlasticity::setUP ( )

private

Definition at line 834 of file multifield_plasticity.cpp.

                                        {
  MoFEMFunctionBegin;
  Simple *simple = mField.getInterface<Simple>();
 
  auto get_base = [this](bool has_hexes = false) {
    MoFEMFunctionBeginHot;
 
    // 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 = has_hexes ? DEMKOWICZ : AINSWORTH;
    if (has_hexes)
      choice_base_value = DEMKOWICZ;
    CHKERR PetscOptionsGetEList(PETSC_NULL, NULL, "-base", list_bases,
                                LASBASETOPT, &choice_base_value, PETSC_NULL);
 
    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;
    }
    MoFEMFunctionReturnHot(0);
  };
 
  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());
  }
 
  CHKERR get_base(nb_hexes);
  // 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) {
    if (data.test_user_base_tau) {
      CHKERR simple->addDomainField("TAU", space_test, USER_BASE, 1);
      CHKERR simple->setFieldOrder("TAU", std::max(0, data.order - 1));
      CHKERR simple->addDomainField("EP", space_test, USER_BASE, 6);
      CHKERR simple->setFieldOrder("EP", std::max(0, data.order - 1));
 
      auto *pipeline_mng = mField.getInterface<PipelineManager>();
 
      // this creates domain_elements
      pipeline_mng->setDomainLhsIntegrationRule(
          [](int, int, int approx_order) { return 1; });
      pipeline_mng->setDomainRhsIntegrationRule(
          [](int, int, int approx_order) { return 1; });
 
      // get access to "TAU" data structure
      for (string name : {"TAU", "EP"}) {
 
        auto field_ptr = mField.get_field_structure(name);
        // get table associating number of dofs to entities depending on
        // approximation order set on those entities.
        auto field_order_table =
            const_cast<Field *>(field_ptr)->getFieldOrderTable();
 
        // function set zero number of dofs
        auto get_tau_bubble_order_zero = [](int p) { return 0; };
        // function set non-zero number of dofs on tetrahedrons
        auto get_tau_bubble_order_tet = [](int p) -> int {
          return std::min(p + 1, 2);
        };
 
        field_order_table[MBVERTEX] = get_tau_bubble_order_zero;
        field_order_table[MBEDGE] = get_tau_bubble_order_zero;
        field_order_table[MBTRI] = get_tau_bubble_order_zero;
        field_order_table[MBTET] = get_tau_bubble_order_tet;
        const_cast<Field *>(field_ptr)->rebuildDofsOrderMap();
      }
 
    } else {
      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 (nb_hexes) {
    CHKERR mField.set_field_order(0, MBQUAD, "U", 0);
    CHKERR mField.set_field_order(0, MBHEX, "U", 0);
    CHKERR mField.build_fields();
  }
  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();
  if (data.with_plasticity)
    CHKERR simple->addFieldToEmptyFieldBlocks("U", "TAU");
 
  MoFEMFunctionReturn(0);
}

tsSolve()

MoFEMErrorCode ContactPlasticity::tsSolve ( )

private

Definition at line 1935 of file multifield_plasticity.cpp.

                                          {
  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, 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)));
    //   if (!roller.directionDataParamName.empty())
    //     preProc->methodsOpForRollersDirection.push_back(
    //         boost::shared_ptr<MethodForForceScaling>(
    //             new TimeAccelerogram(roller.directionDataParamName)));
    // }
 
    // 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_reduced_integration)
        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,
                                     dirichlet_bc_ptr, null);
        CHKERR DMMoFEMTSSetIFunction(dm, simple->getDomainFEName(), null,
                                     preProc, null);
        if (data.is_reduced_integration)
          CHKERR DMMoFEMTSSetIFunction(dm, simple->getDomainFEName(),
                                       feLambdaRhs, null, 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);
 
      auto block_prefix = simple->getProblemName();
      auto bc_mng = mField.getInterface<BcManager>();
 
      auto is_all_bc =
          bc_mng->getBlockIS(block_prefix, "MOVE_Z", "U", name_prb, 2, 2);
      is_all_bc = bc_mng->getBlockIS(block_prefix, "MOVE_Y", "U", name_prb, 1,
                                     1, is_all_bc);
      is_all_bc = bc_mng->getBlockIS(block_prefix, "MOVE_X", "U", name_prb, 0,
                                     0, is_all_bc);
      is_all_bc = bc_mng->getBlockIS(block_prefix, "ROTATE_X", "U", name_prb, 1,
                                     2, is_all_bc);
      is_all_bc = bc_mng->getBlockIS(block_prefix, "ROTATE_Z", "U", name_prb, 0,
                                     1, is_all_bc);
      is_all_bc = bc_mng->getBlockIS(block_prefix, "ROTATE_ALL", "U", name_prb,
                                     0, 2, is_all_bc);
      is_all_bc = bc_mng->getBlockIS(block_prefix, "MOVE_ALL", "U", name_prb, 0,
                                     2, is_all_bc);
 
      int is_all_bc_size;
      CHKERR ISGetSize(is_all_bc, &is_all_bc_size);
 
      if (is_bc_field_split && is_all_bc_size) {
        MOFEM_LOG("WORLD", Sev::inform)
            << "Fieldsplit preconditioner for boundary block size: "
            << is_all_bc_size;
        CHKERR PCFieldSplitSetIS(bc_pc, NULL,
                                 is_all_bc); // boundary block
        //FIXME: it could be beneficial to have 'bc' and 'non_bc' in naming
        
        // CHKERR PCFieldSplitSetIS(bc_pc, "bc",
        //                          is_all_bc); // boundary block
        // IS is_no_bc;
        // CHKERR PCFieldSplitGetISByIndex(bc_pc, 1, &is_no_bc);
        // CHKERR PCFieldSplitGetISByIndex(bc_pc, 0, &is_no_bc);
        // PetscBool is_equal;
        // CHKERR ISEqual(IS is1, IS is2, PetscBool * flg);
        // CHKERR PCFieldSplitSetIS(bc_pc, "non_bc",
        //                          is_no_bc); // non-boundary block
        // CHKERR ISView(is_no_bc, PETSC_VIEWER_STDOUT_SELF);
        // CHKERR ISDestroy(&is_no_bc);
        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 = PETSC_TRUE; // level0 fieldsplit
      PetscBool is_l1_field_split = PETSC_FALSE; // 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 AOView(sigma_ao, PETSC_VIEWER_STDOUT_SELF);
          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);
      };
 
      //FIXME: contact with fieldsplit does not work in parallel, requires further investigation
      //TODO: will come back into this in the future 
      if (data.with_contact && is_l0_field_split && false)
        CHKERR set_fieldsplit_contact();
 
      PetscObjectTypeCompare((PetscObject)pC, PCFIELDSPLIT, &is_l1_field_split);
      if (data.with_plasticity && 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, "ep", is_map["EP"]);
          CHKERR PCFieldSplitSetIS(pC, "etau", 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, "tau", is_map["TAU"]);
              CHKERR PCFieldSplitSetIS(tau_pc, "elastic", 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_20 0 "
                             "-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_20 0 "
                             "-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();
    if (data.test_user_base_tau)
      postProcFe->getUserPolynomialBase() =
          boost::shared_ptr<BaseFunction>(new BubbleTauPolynomialBase());
    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);
}

Member Data Documentation

base

FieldApproximationBase ContactPlasticity::base

Definition at line 157 of file multifield_plasticity.cpp.

bcData

boost::ptr_vector<DataFromMove> ContactPlasticity::bcData

private

Definition at line 141 of file multifield_plasticity.cpp.

boundaryMarker

boost::shared_ptr<std::vector<unsigned char> > ContactPlasticity::boundaryMarker

private

Definition at line 139 of file multifield_plasticity.cpp.

commonDataPtr

boost::shared_ptr<OpContactTools::CommonData> ContactPlasticity::commonDataPtr

private

Definition at line 127 of file multifield_plasticity.cpp.

data

ContactPlasticity::ProblemData ContactPlasticity::data

static

Definition at line 244 of file multifield_plasticity.cpp.

dirichlet_bc_ptr

boost::shared_ptr<DirichletDisplacementRemoveDofsBc> ContactPlasticity::dirichlet_bc_ptr

private

Definition at line 137 of file multifield_plasticity.cpp.

dmContact

SmartPetscObj<DM> ContactPlasticity::dmContact

private

Definition at line 154 of file multifield_plasticity.cpp.

dmEP

SmartPetscObj<DM> ContactPlasticity::dmEP

private

Definition at line 152 of file multifield_plasticity.cpp.

dmTAU

SmartPetscObj<DM> ContactPlasticity::dmTAU

private

Definition at line 153 of file multifield_plasticity.cpp.

edge_forces

boost::ptr_map<std::string, EdgeForce> ContactPlasticity::edge_forces

private

Definition at line 136 of file multifield_plasticity.cpp.

feLambdaLhs

boost::shared_ptr<DomainEle> ContactPlasticity::feLambdaLhs

private

Definition at line 147 of file multifield_plasticity.cpp.

feLambdaRhs

boost::shared_ptr<DomainEle> ContactPlasticity::feLambdaRhs

private

Definition at line 146 of file multifield_plasticity.cpp.

invJac

MatrixDouble ContactPlasticity::invJac

private

Definition at line 126 of file multifield_plasticity.cpp.

jAc

MatrixDouble ContactPlasticity::jAc

private

Definition at line 126 of file multifield_plasticity.cpp.

mField

MoFEM::Interface& ContactPlasticity::mField

private

Definition at line 117 of file multifield_plasticity.cpp.

neumann_forces

boost::ptr_map<std::string, NeumannForcesSurface> ContactPlasticity::neumann_forces

private

Definition at line 134 of file multifield_plasticity.cpp.

nodal_forces

boost::ptr_map<std::string, NodalForce> ContactPlasticity::nodal_forces

private

Definition at line 135 of file multifield_plasticity.cpp.

postProcFe

boost::shared_ptr<PostProcVolumeOnRefinedMesh> ContactPlasticity::postProcFe

private

Definition at line 128 of file multifield_plasticity.cpp.

uXScatter

std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter> > ContactPlasticity::uXScatter

private

Definition at line 129 of file multifield_plasticity.cpp.

uYScatter

std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter> > ContactPlasticity::uYScatter

private

Definition at line 130 of file multifield_plasticity.cpp.

uZScatter

std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter> > ContactPlasticity::uZScatter

private

Definition at line 131 of file multifield_plasticity.cpp.

volSideFe

boost::shared_ptr<DomainSideEle> ContactPlasticity::volSideFe

private

Definition at line 142 of file multifield_plasticity.cpp.

---

The documentation for this struct was generated from the following file:

• users_modules/multifield_plasticity/multifield_plasticity.cpp