843  auto set_section_monitor = [&](
auto solver) {
 
  846    CHKERR TSGetSNES(solver, &snes);
 
  847    CHKERR SNESMonitorSet(snes,
 
  850                          (
void *)(snes_ctx_ptr.get()), 
nullptr);
 
  854  auto create_post_process_elements = [&]() {
 
  855    auto push_vol_ops = [
this](
auto &pip) {
 
  857          pip, {
H1, 
HDIV}, 
"GEOMETRY");
 
  859      auto [common_plastic_ptr, common_hencky_ptr] =
 
  860          PlasticOps::createCommonPlasticOps<SPACE_DIM, IT, DomainEleOp>(
 
  861              mField, 
"MAT_PLASTIC", pip, 
"U", 
"EP", 
"TAU", 1., Sev::inform);
 
  863      if (common_hencky_ptr) {
 
  864        if (common_plastic_ptr->mGradPtr != common_hencky_ptr->matGradPtr)
 
  868      return std::make_pair(common_plastic_ptr, common_hencky_ptr);
 
  871    auto push_vol_post_proc_ops = [
this](
auto &pp_fe, 
auto &&p) {
 
  874      auto &pip = pp_fe->getOpPtrVector();
 
  876      auto [common_plastic_ptr, common_hencky_ptr] = p;
 
  880      auto x_ptr = boost::make_shared<MatrixDouble>();
 
  883      auto u_ptr = boost::make_shared<MatrixDouble>();
 
  892                pp_fe->getPostProcMesh(), pp_fe->getMapGaussPts(),
 
  895                  common_plastic_ptr->getPlasticSurfacePtr()},
 
  896                 {
"PLASTIC_MULTIPLIER",
 
  897                  common_plastic_ptr->getPlasticTauPtr()}},
 
  899                {{
"U", u_ptr}, {
"GEOMETRY", x_ptr}},
 
  901                {{
"GRAD", common_hencky_ptr->matGradPtr},
 
  902                 {
"FIRST_PIOLA", common_hencky_ptr->getMatFirstPiolaStress()}},
 
  904                {{
"HENCKY_STRAIN", common_hencky_ptr->getMatLogC()},
 
  905                 {
"PLASTIC_STRAIN", common_plastic_ptr->getPlasticStrainPtr()},
 
  906                 {
"PLASTIC_FLOW", common_plastic_ptr->getPlasticFlowPtr()}}
 
  918                pp_fe->getPostProcMesh(), pp_fe->getMapGaussPts(),
 
  921                  common_plastic_ptr->getPlasticSurfacePtr()},
 
  922                 {
"PLASTIC_MULTIPLIER",
 
  923                  common_plastic_ptr->getPlasticTauPtr()}},
 
  925                {{
"U", u_ptr}, {
"GEOMETRY", x_ptr}},
 
  929                {{
"STRAIN", common_plastic_ptr->mStrainPtr},
 
  930                 {
"STRESS", common_plastic_ptr->mStressPtr},
 
  931                 {
"PLASTIC_STRAIN", common_plastic_ptr->getPlasticStrainPtr()},
 
  932                 {
"PLASTIC_FLOW", common_plastic_ptr->getPlasticFlowPtr()}}
 
  942    PetscBool post_proc_vol;
 
  943    PetscBool post_proc_skin;
 
  946      post_proc_vol = PETSC_TRUE;
 
  947      post_proc_skin = PETSC_FALSE;
 
  949      post_proc_vol = PETSC_FALSE;
 
  950      post_proc_skin = PETSC_TRUE;
 
  955                               &post_proc_skin, PETSC_NULLPTR);
 
  957    auto vol_post_proc = [
this, push_vol_post_proc_ops, push_vol_ops,
 
  959      if (post_proc_vol == PETSC_FALSE)
 
  960        return boost::shared_ptr<PostProcEle>();
 
  961      auto pp_fe = boost::make_shared<PostProcEle>(mField);
 
  963          push_vol_post_proc_ops(pp_fe, push_vol_ops(pp_fe->getOpPtrVector())),
 
  964          "push_vol_post_proc_ops");
 
  968    auto skin_post_proc = [
this, push_vol_post_proc_ops, push_vol_ops,
 
  970      if (post_proc_skin == PETSC_FALSE)
 
  971        return boost::shared_ptr<SkinPostProcEle>();
 
  974      auto pp_fe = boost::make_shared<SkinPostProcEle>(mField);
 
  977      pp_fe->getOpPtrVector().push_back(op_side);
 
  979                         pp_fe, push_vol_ops(op_side->getOpPtrVector())),
 
  980                     "push_vol_post_proc_ops");
 
  984    return std::make_pair(vol_post_proc(), skin_post_proc());
 
  987  auto scatter_create = [&](
auto D, 
auto coeff) {
 
  989    CHKERR is_manager->isCreateProblemFieldAndRank(
simple->getProblemName(),
 
  990                                                   ROW, 
"U", coeff, coeff, is);
 
  992    CHKERR ISGetLocalSize(is, &loc_size);
 
  994    CHKERR VecCreateMPI(mField.get_comm(), loc_size, PETSC_DETERMINE, &
v);
 
  996    CHKERR VecScatterCreate(
D, is, 
v, PETSC_NULLPTR, &scatter);
 
 1001  boost::shared_ptr<SetPtsData> field_eval_data;
 
 1002  boost::shared_ptr<MatrixDouble> u_field_ptr;
 
 1004  std::array<double, 3> field_eval_coords{0.0, 0.0, 0.0};
 
 1007                                  field_eval_coords.data(), &coords_dim,
 
 1010  boost::shared_ptr<std::map<std::string, boost::shared_ptr<VectorDouble>>>
 
 1011      scalar_field_ptrs = boost::make_shared<
 
 1012          std::map<std::string, boost::shared_ptr<VectorDouble>>>();
 
 1013  boost::shared_ptr<std::map<std::string, boost::shared_ptr<MatrixDouble>>>
 
 1014      vector_field_ptrs = boost::make_shared<
 
 1015          std::map<std::string, boost::shared_ptr<MatrixDouble>>>();
 
 1016  boost::shared_ptr<std::map<std::string, boost::shared_ptr<MatrixDouble>>>
 
 1017      sym_tensor_field_ptrs = boost::make_shared<
 
 1018          std::map<std::string, boost::shared_ptr<MatrixDouble>>>();
 
 1019  boost::shared_ptr<std::map<std::string, boost::shared_ptr<MatrixDouble>>>
 
 1020      tensor_field_ptrs = boost::make_shared<
 
 1021          std::map<std::string, boost::shared_ptr<MatrixDouble>>>();
 
 1024    auto u_field_ptr = boost::make_shared<MatrixDouble>();
 
 1029        field_eval_data, 
simple->getDomainFEName());
 
 1031    field_eval_data->setEvalPoints(field_eval_coords.data(), 1);
 
 1032    auto no_rule = [](
int, 
int, 
int) { 
return -1; };
 
 1033    auto field_eval_fe_ptr = field_eval_data->feMethodPtr.lock();
 
 1034    field_eval_fe_ptr->getRuleHook = no_rule;
 
 1037        field_eval_fe_ptr->getOpPtrVector(), {H1, HDIV}, 
"GEOMETRY");
 
 1039    auto [common_plastic_ptr, common_hencky_ptr] =
 
 1040        PlasticOps::createCommonPlasticOps<SPACE_DIM, IT, DomainEleOp>(
 
 1041            mField, 
"MAT_PLASTIC", field_eval_fe_ptr->getOpPtrVector(), 
"U",
 
 1042            "EP", 
"TAU", 1., Sev::inform);
 
 1044    field_eval_fe_ptr->getOpPtrVector().push_back(
 
 1047    if ((common_plastic_ptr) && (common_hencky_ptr) && (scalar_field_ptrs)) {
 
 1049        scalar_field_ptrs->insert(
 
 1050            {
"PLASTIC_SURFACE", common_plastic_ptr->getPlasticSurfacePtr()});
 
 1051        scalar_field_ptrs->insert(
 
 1052            {
"PLASTIC_MULTIPLIER", common_plastic_ptr->getPlasticTauPtr()});
 
 1053        vector_field_ptrs->insert({
"U", u_field_ptr});
 
 1054        sym_tensor_field_ptrs->insert(
 
 1055            {
"PLASTIC_STRAIN", common_plastic_ptr->getPlasticStrainPtr()});
 
 1056        sym_tensor_field_ptrs->insert(
 
 1057            {
"PLASTIC_FLOW", common_plastic_ptr->getPlasticFlowPtr()});
 
 1058        sym_tensor_field_ptrs->insert(
 
 1059            {
"HENCKY_STRAIN", common_hencky_ptr->getMatLogC()});
 
 1060        tensor_field_ptrs->insert({
"GRAD", common_hencky_ptr->matGradPtr});
 
 1061        tensor_field_ptrs->insert(
 
 1062            {
"FIRST_PIOLA", common_hencky_ptr->getMatFirstPiolaStress()});
 
 1064        scalar_field_ptrs->insert(
 
 1065            {
"PLASTIC_SURFACE", common_plastic_ptr->getPlasticSurfacePtr()});
 
 1066        scalar_field_ptrs->insert(
 
 1067            {
"PLASTIC_MULTIPLIER", common_plastic_ptr->getPlasticTauPtr()});
 
 1068        vector_field_ptrs->insert({
"U", u_field_ptr});
 
 1069        sym_tensor_field_ptrs->insert(
 
 1070            {
"STRAIN", common_plastic_ptr->mStrainPtr});
 
 1071        sym_tensor_field_ptrs->insert(
 
 1072            {
"STRESS", common_plastic_ptr->mStressPtr});
 
 1073        sym_tensor_field_ptrs->insert(
 
 1074            {
"PLASTIC_STRAIN", common_plastic_ptr->getPlasticStrainPtr()});
 
 1075        sym_tensor_field_ptrs->insert(
 
 1076            {
"PLASTIC_FLOW", common_plastic_ptr->getPlasticFlowPtr()});
 
 1081  auto test_monitor_ptr = boost::make_shared<FEMethod>();
 
 1083  auto set_time_monitor = [&](
auto dm, 
auto solver) {
 
 1086        dm, create_post_process_elements(), reactionFe, uXScatter, uYScatter,
 
 1087        uZScatter, field_eval_coords, field_eval_data, scalar_field_ptrs,
 
 1088        vector_field_ptrs, sym_tensor_field_ptrs, tensor_field_ptrs));
 
 1089    boost::shared_ptr<ForcesAndSourcesCore> null;
 
 1091    test_monitor_ptr->postProcessHook = [&]() {
 
 1094      if (
atom_test && fabs(test_monitor_ptr->ts_t - 0.5) < 1e-12 &&
 
 1095          test_monitor_ptr->ts_step == 25) {
 
 1097        if (scalar_field_ptrs->at(
"PLASTIC_MULTIPLIER")->size()) {
 
 1100          MOFEM_LOG(
"PlasticSync", Sev::inform) << 
"Eval point tau: " << t_tau;
 
 1102          if (
atom_test == 1 && fabs(t_tau - 0.688861) > 1e-5) {
 
 1104                     "atom test %d failed: wrong plastic multiplier value",
 
 1109        if (vector_field_ptrs->at(
"U")->size1()) {
 
 1112              getFTensor1FromMat<SPACE_DIM>(*vector_field_ptrs->at(
"U"));
 
 1113          MOFEM_LOG(
"PlasticSync", Sev::inform) << 
"Eval point U: " << t_disp;
 
 1115          if (
atom_test == 1 && fabs(t_disp(0) - 0.25 / 2.) > 1e-5 ||
 
 1116              fabs(t_disp(1) + 0.0526736) > 1e-5) {
 
 1118                     "atom test %d failed: wrong displacement value",
 
 1123        if (sym_tensor_field_ptrs->at(
"PLASTIC_STRAIN")->size1()) {
 
 1124          auto t_plastic_strain = getFTensor2SymmetricFromMat<SPACE_DIM>(
 
 1125              *sym_tensor_field_ptrs->at(
"PLASTIC_STRAIN"));
 
 1127              << 
"Eval point EP: " << t_plastic_strain;
 
 1130                  fabs(t_plastic_strain(0, 0) - 0.221943) > 1e-5 ||
 
 1131              fabs(t_plastic_strain(0, 1)) > 1e-5 ||
 
 1132              fabs(t_plastic_strain(1, 1) + 0.110971) > 1e-5) {
 
 1134                     "atom test %d failed: wrong plastic strain value",
 
 1139        if (tensor_field_ptrs->at(
"FIRST_PIOLA")->size1()) {
 
 1140          auto t_piola_stress = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(
 
 1141              *tensor_field_ptrs->at(
"FIRST_PIOLA"));
 
 1143              << 
"Eval point Piola stress: " << t_piola_stress;
 
 1145          if (
atom_test == 1 && fabs((t_piola_stress(0, 0) - 198.775) /
 
 1146                                     t_piola_stress(0, 0)) > 1e-5 ||
 
 1147              fabs(t_piola_stress(0, 1)) + fabs(t_piola_stress(1, 0)) +
 
 1148                      fabs(t_piola_stress(1, 1)) >
 
 1151                     "atom test %d failed: wrong Piola stress value",
 
 1162                               monitor_ptr, null, test_monitor_ptr);
 
 1167  auto set_schur_pc = [&](
auto solver,
 
 1168                          boost::shared_ptr<SetUpSchur> &schur_ptr) {
 
 1171    auto name_prb = 
simple->getProblemName();
 
 1177      dm_sub = 
createDM(mField.get_comm(), 
"DMMOFEM");
 
 1182      for (
auto f : {
"U"}) {
 
 1194      dm_sub = 
createDM(mField.get_comm(), 
"DMMOFEM");
 
 1200      for (
auto f : {
"SIGMA", 
"EP", 
"TAU"}) {
 
 1205      for (
auto f : {
"EP", 
"TAU"}) {
 
 1215    if constexpr (
AT == AssemblyType::BLOCK_SCHUR) {
 
 1224      auto get_nested_mat_data = [&](
auto schur_dm, 
auto block_dm) {
 
 1230                {simple->getDomainFEName(),
 
 1246                {simple->getBoundaryFEName(),
 
 1248                 {{
"SIGMA", 
"SIGMA"}, {
"U", 
"SIGMA"}, {
"SIGMA", 
"U"}
 
 1258            {dm_schur, dm_block}, block_mat_data,
 
 1260            {
"SIGMA", 
"EP", 
"TAU"}, {
nullptr, 
nullptr, 
nullptr}, true
 
 1267      auto get_nested_mat_data = [&](
auto schur_dm, 
auto block_dm) {
 
 1268        auto block_mat_data =
 
 1271                                    {{simple->getDomainFEName(),
 
 1288            {dm_schur, dm_block}, block_mat_data,
 
 1290            {
"EP", 
"TAU"}, {
nullptr, 
nullptr}, false
 
 1297      auto nested_mat_data = get_nested_mat_data(dm_schur, dm_block);
 
 1300      auto block_is = 
getDMSubData(dm_block)->getSmartRowIs();
 
 1301      auto ao_schur = 
getDMSubData(dm_schur)->getSmartRowMap();
 
 1307      CHKERR schur_ptr->setUp(solver);
 
 1313  auto dm = 
simple->getDM();
 
 1316  uXScatter = scatter_create(
D, 0);
 
 1317  uYScatter = scatter_create(
D, 1);
 
 1319    uZScatter = scatter_create(
D, 2);
 
 1321  auto create_solver = [pip_mng]() {
 
 1323      return pip_mng->createTSIM();
 
 1325      return pip_mng->createTSIM2();
 
 1328  auto solver = create_solver();
 
 1330  auto active_pre_lhs = []() {
 
 1337  auto active_post_lhs = [&]() {
 
 1339    auto get_iter = [&]() {
 
 1344                        "Can not get iter");
 
 1348    auto iter = get_iter();
 
 1351      std::array<int, 5> activity_data;
 
 1352      std::fill(activity_data.begin(), activity_data.end(), 0);
 
 1354                    activity_data.data(), activity_data.size(), MPI_INT,
 
 1355                    MPI_SUM, mField.get_comm());
 
 1357      int &active_points = activity_data[0];
 
 1358      int &avtive_full_elems = activity_data[1];
 
 1359      int &avtive_elems = activity_data[2];
 
 1360      int &nb_points = activity_data[3];
 
 1361      int &nb_elements = activity_data[4];
 
 1365        double proc_nb_points =
 
 1366            100 * 
static_cast<double>(active_points) / nb_points;
 
 1367        double proc_nb_active =
 
 1368            100 * 
static_cast<double>(avtive_elems) / nb_elements;
 
 1369        double proc_nb_full_active = 100;
 
 1371          proc_nb_full_active =
 
 1372              100 * 
static_cast<double>(avtive_full_elems) / avtive_elems;
 
 1375                    "Iter %d nb pts %d nb active pts %d (%3.3f\%) nb active " 
 1377                    "(%3.3f\%) nb full active elems %d (%3.3f\%)",
 
 1378                    iter, nb_points, active_points, proc_nb_points,
 
 1379                    avtive_elems, proc_nb_active, avtive_full_elems,
 
 1380                    proc_nb_full_active, iter);
 
 1387  auto add_active_dofs_elem = [&](
auto dm) {
 
 1389    auto fe_pre_proc = boost::make_shared<FEMethod>();
 
 1390    fe_pre_proc->preProcessHook = active_pre_lhs;
 
 1391    auto fe_post_proc = boost::make_shared<FEMethod>();
 
 1392    fe_post_proc->postProcessHook = active_post_lhs;
 
 1394    ts_ctx_ptr->getPreProcessIJacobian().push_front(fe_pre_proc);
 
 1395    ts_ctx_ptr->getPostProcessIJacobian().push_back(fe_post_proc);
 
 1399  auto set_essential_bc = [&](
auto dm, 
auto solver) {
 
 1403    auto pre_proc_ptr = boost::make_shared<FEMethod>();
 
 1404    auto post_proc_rhs_ptr = boost::make_shared<FEMethod>();
 
 1405    auto post_proc_lhs_ptr = boost::make_shared<FEMethod>();
 
 1408    auto disp_time_scale = boost::make_shared<TimeScale>();
 
 1410    auto get_bc_hook_rhs = [
this, pre_proc_ptr, disp_time_scale]() {
 
 1412                                                     {disp_time_scale}, 
false);
 
 1415    pre_proc_ptr->preProcessHook = get_bc_hook_rhs();
 
 1417    auto get_post_proc_hook_rhs = [
this, post_proc_rhs_ptr]() {
 
 1420          mField, post_proc_rhs_ptr, 
nullptr, Sev::verbose)();
 
 1422          mField, post_proc_rhs_ptr, 1.)();
 
 1425    auto get_post_proc_hook_lhs = [
this, post_proc_lhs_ptr]() {
 
 1427          mField, post_proc_lhs_ptr, 1.);
 
 1429    post_proc_rhs_ptr->postProcessHook = get_post_proc_hook_rhs;
 
 1432    ts_ctx_ptr->getPreProcessIFunction().push_front(pre_proc_ptr);
 
 1433    ts_ctx_ptr->getPreProcessIJacobian().push_front(pre_proc_ptr);
 
 1434    ts_ctx_ptr->getPostProcessIFunction().push_back(post_proc_rhs_ptr);
 
 1435    post_proc_lhs_ptr->postProcessHook = get_post_proc_hook_lhs();
 
 1436    ts_ctx_ptr->getPostProcessIJacobian().push_back(post_proc_lhs_ptr);
 
 1443    CHKERR TSSetIJacobian(solver, 
B, 
B, PETSC_NULLPTR, PETSC_NULLPTR);
 
 1445    CHKERR TSSetI2Jacobian(solver, 
B, 
B, PETSC_NULLPTR, PETSC_NULLPTR);
 
 1448    CHKERR TSSetSolution(solver, 
D);
 
 1450    CHKERR TS2SetSolution(solver, 
D, DD);
 
 1452  CHKERR set_section_monitor(solver);
 
 1453  CHKERR set_time_monitor(dm, solver);
 
 1454  CHKERR TSSetFromOptions(solver);
 
 1456  CHKERR add_active_dofs_elem(dm);
 
 1457  boost::shared_ptr<SetUpSchur> schur_ptr;
 
 1458  CHKERR set_schur_pc(solver, schur_ptr);
 
 1459  CHKERR set_essential_bc(dm, solver);
 
 1464    BOOST_LOG_SCOPED_THREAD_ATTR(
"Timeline", attrs::timer());
 
 1465  MOFEM_LOG(
"TIMER", Sev::verbose) << 
"TSSetUp";
 
 1467  MOFEM_LOG(
"TIMER", Sev::verbose) << 
"TSSetUp <= done";
 
 1468  MOFEM_LOG(
"TIMER", Sev::verbose) << 
"TSSolve";
 
 1469  CHKERR TSSolve(solver, NULL);
 
 1470  MOFEM_LOG(
"TIMER", Sev::verbose) << 
"TSSolve <= done";