835 {
837
841
843
844 auto set_section_monitor = [&](auto solver) {
846 SNES snes;
847 CHKERR TSGetSNES(solver, &snes);
848 CHKERR SNESMonitorSet(snes,
851 (void *)(snes_ctx_ptr.get()), nullptr);
853 };
854
855 auto create_post_process_elements = [&]() {
856 auto push_vol_ops = [this](auto &pip) {
858 pip, {
H1,
HDIV},
"GEOMETRY");
859
860 auto [common_plastic_ptr, common_hencky_ptr] =
861 PlasticOps::createCommonPlasticOps<SPACE_DIM, IT, DomainEleOp>(
862 mField,
"MAT_PLASTIC", pip,
"U",
"EP",
"TAU", 1., Sev::inform);
863
864 if (common_hencky_ptr) {
865 if (common_plastic_ptr->mGradPtr != common_hencky_ptr->matGradPtr)
867 }
868
869 return std::make_pair(common_plastic_ptr, common_hencky_ptr);
870 };
871
872 auto push_vol_post_proc_ops = [this](auto &pp_fe, auto &&p) {
874
875 auto &pip = pp_fe->getOpPtrVector();
876
877 auto [common_plastic_ptr, common_hencky_ptr] = p;
878
880
881 auto x_ptr = boost::make_shared<MatrixDouble>();
882 pip.push_back(
884 auto u_ptr = boost::make_shared<MatrixDouble>();
886
888
889 pip.push_back(
890
892
893 pp_fe->getPostProcMesh(), pp_fe->getMapGaussPts(),
894
895 {{"PLASTIC_SURFACE",
896 common_plastic_ptr->getPlasticSurfacePtr()},
897 {"PLASTIC_MULTIPLIER",
898 common_plastic_ptr->getPlasticTauPtr()}},
899
900 {{"U", u_ptr}, {"GEOMETRY", x_ptr}},
901
902 {{"GRAD", common_hencky_ptr->matGradPtr},
903 {"FIRST_PIOLA", common_hencky_ptr->getMatFirstPiolaStress()}},
904
905 {{"HENCKY_STRAIN", common_hencky_ptr->getMatLogC()},
906 {"PLASTIC_STRAIN", common_plastic_ptr->getPlasticStrainPtr()},
907 {"PLASTIC_FLOW", common_plastic_ptr->getPlasticFlowPtr()}}
908
909 )
910
911 );
912
913 } else {
914
915 pip.push_back(
916
918
919 pp_fe->getPostProcMesh(), pp_fe->getMapGaussPts(),
920
921 {{"PLASTIC_SURFACE",
922 common_plastic_ptr->getPlasticSurfacePtr()},
923 {"PLASTIC_MULTIPLIER",
924 common_plastic_ptr->getPlasticTauPtr()}},
925
926 {{"U", u_ptr}, {"GEOMETRY", x_ptr}},
927
928 {},
929
930 {{"STRAIN", common_plastic_ptr->mStrainPtr},
931 {"STRESS", common_plastic_ptr->mStressPtr},
932 {"PLASTIC_STRAIN", common_plastic_ptr->getPlasticStrainPtr()},
933 {"PLASTIC_FLOW", common_plastic_ptr->getPlasticFlowPtr()}}
934
935 )
936
937 );
938 }
939
941 };
942
943 PetscBool post_proc_vol;
944 PetscBool post_proc_skin;
945
947 post_proc_vol = PETSC_TRUE;
948 post_proc_skin = PETSC_FALSE;
949 } else {
950 post_proc_vol = PETSC_FALSE;
951 post_proc_skin = PETSC_TRUE;
952 }
954 PETSC_NULLPTR);
956 &post_proc_skin, PETSC_NULLPTR);
957
958 auto vol_post_proc = [this, push_vol_post_proc_ops, push_vol_ops,
959 post_proc_vol]() {
960 if (post_proc_vol == PETSC_FALSE)
961 return boost::shared_ptr<PostProcEle>();
962 auto pp_fe = boost::make_shared<PostProcEle>(
mField);
964 push_vol_post_proc_ops(pp_fe, push_vol_ops(pp_fe->getOpPtrVector())),
965 "push_vol_post_proc_ops");
966 return pp_fe;
967 };
968
969 auto skin_post_proc = [this, push_vol_post_proc_ops, push_vol_ops,
970 post_proc_skin]() {
971 if (post_proc_skin == PETSC_FALSE)
972 return boost::shared_ptr<SkinPostProcEle>();
973
975 auto pp_fe = boost::make_shared<SkinPostProcEle>(
mField);
978 pp_fe->getOpPtrVector().push_back(op_side);
980 pp_fe, push_vol_ops(op_side->getOpPtrVector())),
981 "push_vol_post_proc_ops");
982 return pp_fe;
983 };
984
985 return std::make_pair(vol_post_proc(), skin_post_proc());
986 };
987
988 auto scatter_create = [&](
auto D,
auto coeff) {
990 CHKERR is_manager->isCreateProblemFieldAndRank(
simple->getProblemName(),
991 ROW,
"U", coeff, coeff, is);
992 int loc_size;
993 CHKERR ISGetLocalSize(is, &loc_size);
996 VecScatter scatter;
997 CHKERR VecScatterCreate(
D, is,
v, PETSC_NULLPTR, &scatter);
1000 };
1001
1002 boost::shared_ptr<SetPtsData> field_eval_data;
1003 boost::shared_ptr<MatrixDouble> u_field_ptr;
1004
1005 std::array<double, 3> field_eval_coords{0.0, 0.0, 0.0};
1006 int coords_dim = 3;
1008 field_eval_coords.data(), &coords_dim,
1010
1011 boost::shared_ptr<std::map<std::string, boost::shared_ptr<VectorDouble>>>
1012 scalar_field_ptrs = boost::make_shared<
1013 std::map<std::string, boost::shared_ptr<VectorDouble>>>();
1014 boost::shared_ptr<std::map<std::string, boost::shared_ptr<MatrixDouble>>>
1015 vector_field_ptrs = boost::make_shared<
1016 std::map<std::string, boost::shared_ptr<MatrixDouble>>>();
1017 boost::shared_ptr<std::map<std::string, boost::shared_ptr<MatrixDouble>>>
1018 sym_tensor_field_ptrs = boost::make_shared<
1019 std::map<std::string, boost::shared_ptr<MatrixDouble>>>();
1020 boost::shared_ptr<std::map<std::string, boost::shared_ptr<MatrixDouble>>>
1021 tensor_field_ptrs = boost::make_shared<
1022 std::map<std::string, boost::shared_ptr<MatrixDouble>>>();
1023
1025 auto u_field_ptr = boost::make_shared<MatrixDouble>();
1026 field_eval_data =
1028
1030 field_eval_data,
simple->getDomainFEName());
1031
1032 field_eval_data->setEvalPoints(field_eval_coords.data(), 1);
1033 auto no_rule = [](
int,
int,
int) {
return -1; };
1034 auto field_eval_fe_ptr = field_eval_data->feMethodPtr;
1035 field_eval_fe_ptr->getRuleHook = no_rule;
1036
1038 field_eval_fe_ptr->getOpPtrVector(), {H1, HDIV}, "GEOMETRY");
1039
1040 auto [common_plastic_ptr, common_hencky_ptr] =
1041 PlasticOps::createCommonPlasticOps<SPACE_DIM, IT, DomainEleOp>(
1042 mField,
"MAT_PLASTIC", field_eval_fe_ptr->getOpPtrVector(),
"U",
1043 "EP", "TAU", 1., Sev::inform);
1044
1045 field_eval_fe_ptr->getOpPtrVector().push_back(
1047
1048 if ((common_plastic_ptr) && (common_hencky_ptr) && (scalar_field_ptrs)) {
1050 scalar_field_ptrs->insert(
1051 {"PLASTIC_SURFACE", common_plastic_ptr->getPlasticSurfacePtr()});
1052 scalar_field_ptrs->insert(
1053 {"PLASTIC_MULTIPLIER", common_plastic_ptr->getPlasticTauPtr()});
1054 vector_field_ptrs->insert({"U", u_field_ptr});
1055 sym_tensor_field_ptrs->insert(
1056 {"PLASTIC_STRAIN", common_plastic_ptr->getPlasticStrainPtr()});
1057 sym_tensor_field_ptrs->insert(
1058 {"PLASTIC_FLOW", common_plastic_ptr->getPlasticFlowPtr()});
1059 sym_tensor_field_ptrs->insert(
1060 {"HENCKY_STRAIN", common_hencky_ptr->getMatLogC()});
1061 tensor_field_ptrs->insert({"GRAD", common_hencky_ptr->matGradPtr});
1062 tensor_field_ptrs->insert(
1063 {"FIRST_PIOLA", common_hencky_ptr->getMatFirstPiolaStress()});
1064 } else {
1065 scalar_field_ptrs->insert(
1066 {"PLASTIC_SURFACE", common_plastic_ptr->getPlasticSurfacePtr()});
1067 scalar_field_ptrs->insert(
1068 {"PLASTIC_MULTIPLIER", common_plastic_ptr->getPlasticTauPtr()});
1069 vector_field_ptrs->insert({"U", u_field_ptr});
1070 sym_tensor_field_ptrs->insert(
1071 {"STRAIN", common_plastic_ptr->mStrainPtr});
1072 sym_tensor_field_ptrs->insert(
1073 {"STRESS", common_plastic_ptr->mStressPtr});
1074 sym_tensor_field_ptrs->insert(
1075 {"PLASTIC_STRAIN", common_plastic_ptr->getPlasticStrainPtr()});
1076 sym_tensor_field_ptrs->insert(
1077 {"PLASTIC_FLOW", common_plastic_ptr->getPlasticFlowPtr()});
1078 }
1079 }
1080 }
1081
1082 auto test_monitor_ptr = boost::make_shared<FEMethod>();
1083
1084 auto set_time_monitor = [&](auto dm, auto solver) {
1088 uZScatter, field_eval_coords, field_eval_data, scalar_field_ptrs,
1089 vector_field_ptrs, sym_tensor_field_ptrs, tensor_field_ptrs));
1090 boost::shared_ptr<ForcesAndSourcesCore> null;
1091
1092 test_monitor_ptr->postProcessHook = [&]() {
1094
1095 if (
atom_test && fabs(test_monitor_ptr->ts_t - 0.5) < 1e-12 &&
1096 test_monitor_ptr->ts_step == 25) {
1097
1098 if (scalar_field_ptrs->at("PLASTIC_MULTIPLIER")->size()) {
1099 auto t_tau =
1101 MOFEM_LOG(
"PlasticSync", Sev::inform) <<
"Eval point tau: " << t_tau;
1102
1103 if (
atom_test == 1 && fabs(t_tau - 0.688861) > 1e-5) {
1105 "atom test %d failed: wrong plastic multiplier value",
1107 }
1108 }
1109
1110 if (vector_field_ptrs->at("U")->size1()) {
1112 auto t_disp =
1113 getFTensor1FromMat<SPACE_DIM>(*vector_field_ptrs->at("U"));
1114 MOFEM_LOG(
"PlasticSync", Sev::inform) <<
"Eval point U: " << t_disp;
1115
1116 if (
atom_test == 1 && fabs(t_disp(0) - 0.25 / 2.) > 1e-5 ||
1117 fabs(t_disp(1) + 0.0526736) > 1e-5) {
1119 "atom test %d failed: wrong displacement value",
1121 }
1122 }
1123
1124 if (sym_tensor_field_ptrs->at("PLASTIC_STRAIN")->size1()) {
1125 auto t_plastic_strain = getFTensor2SymmetricFromMat<SPACE_DIM>(
1126 *sym_tensor_field_ptrs->at("PLASTIC_STRAIN"));
1128 << "Eval point EP: " << t_plastic_strain;
1129
1131 fabs(t_plastic_strain(0, 0) - 0.221943) > 1e-5 ||
1132 fabs(t_plastic_strain(0, 1)) > 1e-5 ||
1133 fabs(t_plastic_strain(1, 1) + 0.110971) > 1e-5) {
1135 "atom test %d failed: wrong plastic strain value",
1137 }
1138 }
1139
1140 if (tensor_field_ptrs->at("FIRST_PIOLA")->size1()) {
1141 auto t_piola_stress = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(
1142 *tensor_field_ptrs->at("FIRST_PIOLA"));
1144 << "Eval point Piola stress: " << t_piola_stress;
1145
1146 if (
atom_test == 1 && fabs((t_piola_stress(0, 0) - 198.775) /
1147 t_piola_stress(0, 0)) > 1e-5 ||
1148 fabs(t_piola_stress(0, 1)) + fabs(t_piola_stress(1, 0)) +
1149 fabs(t_piola_stress(1, 1)) >
1150 1e-5) {
1152 "atom test %d failed: wrong Piola stress value",
1154 }
1155 }
1156 }
1157
1160 };
1161
1163 monitor_ptr, null, test_monitor_ptr);
1164
1166 };
1167
1168 auto set_schur_pc = [&](auto solver,
1169 boost::shared_ptr<SetUpSchur> &schur_ptr) {
1171
1172 auto name_prb =
simple->getProblemName();
1173
1174
1183 for (
auto f : {
"U"}) {
1186 }
1188
1190 };
1191
1200#ifdef ADD_CONTACT
1201 for (
auto f : {
"SIGMA",
"EP",
"TAU"}) {
1204 }
1205#else
1206 for (
auto f : {
"EP",
"TAU"}) {
1209 }
1210#endif
1213 };
1214
1215
1216 if constexpr (
AT == AssemblyType::BLOCK_SCHUR) {
1217
1222
1223#ifdef ADD_CONTACT
1224
1225 auto get_nested_mat_data = [&](auto schur_dm, auto block_dm) {
1228
1229 {
1230
1231 {simple->getDomainFEName(),
1232
1233 {{"U", "U"},
1234 {"SIGMA", "SIGMA"},
1235 {"U", "SIGMA"},
1236 {"SIGMA", "U"},
1237 {"EP", "EP"},
1238 {"TAU", "TAU"},
1239 {"U", "EP"},
1240 {"EP", "U"},
1241 {"EP", "TAU"},
1242 {"TAU", "EP"},
1243 {"TAU", "U"}
1244
1245 }},
1246
1247 {simple->getBoundaryFEName(),
1248
1249 {{"SIGMA", "SIGMA"}, {"U", "SIGMA"}, {"SIGMA", "U"}
1250
1251 }}
1252
1253 }
1254
1255 );
1256
1258
1259 {dm_schur, dm_block}, block_mat_data,
1260
1261 {"SIGMA", "EP", "TAU"}, {nullptr, nullptr, nullptr}, true
1262
1263 );
1264 };
1265
1266#else
1267
1268 auto get_nested_mat_data = [&](auto schur_dm, auto block_dm) {
1269 auto block_mat_data =
1271
1272 {{simple->getDomainFEName(),
1273
1274 {{"U", "U"},
1275 {"EP", "EP"},
1276 {"TAU", "TAU"},
1277 {"U", "EP"},
1278 {"EP", "U"},
1279 {"EP", "TAU"},
1280 {"TAU", "U"},
1281 {"TAU", "EP"}
1282
1283 }}}
1284
1285 );
1286
1288
1289 {dm_schur, dm_block}, block_mat_data,
1290
1291 {"EP", "TAU"}, {nullptr, nullptr}, false
1292
1293 );
1294 };
1295
1296#endif
1297
1298 auto nested_mat_data = get_nested_mat_data(dm_schur, dm_block);
1300
1301 auto block_is =
getDMSubData(dm_block)->getSmartRowIs();
1302 auto ao_schur =
getDMSubData(dm_schur)->getSmartRowMap();
1303
1304
1305
1306 schur_ptr =
1308 CHKERR schur_ptr->setUp(solver);
1309 }
1310
1312 };
1313
1314 auto dm =
simple->getDM();
1317 CHKERR VecSetDM(
D, PETSC_NULLPTR);
1318 CHKERR VecSetDM(DD, PETSC_NULLPTR);
1323
1324 auto create_solver = [pip_mng]() {
1326 return pip_mng->createTSIM();
1327 else
1328 return pip_mng->createTSIM2();
1329 };
1330
1331 auto solver = create_solver();
1332
1333 auto active_pre_lhs = []() {
1338 };
1339
1340 auto active_post_lhs = [&]() {
1342 auto get_iter = [&]() {
1343 SNES snes;
1345 int iter;
1347 "Can not get iter");
1348 return iter;
1349 };
1350
1351 auto iter = get_iter();
1352 if (iter >= 0) {
1353
1354 std::array<int, 5> activity_data;
1355 std::fill(activity_data.begin(), activity_data.end(), 0);
1357 activity_data.data(), activity_data.size(), MPI_INT,
1359
1360 int &active_points = activity_data[0];
1361 int &avtive_full_elems = activity_data[1];
1362 int &avtive_elems = activity_data[2];
1363 int &nb_points = activity_data[3];
1364 int &nb_elements = activity_data[4];
1365
1366 if (nb_points) {
1367
1368 double proc_nb_points =
1369 100 * static_cast<double>(active_points) / nb_points;
1370 double proc_nb_active =
1371 100 * static_cast<double>(avtive_elems) / nb_elements;
1372 double proc_nb_full_active = 100;
1373 if (avtive_elems)
1374 proc_nb_full_active =
1375 100 * static_cast<double>(avtive_full_elems) / avtive_elems;
1376
1378 "Iter %d nb pts %d nb active pts %d (%3.3f\%) nb active "
1379 "elements %d "
1380 "(%3.3f\%) nb full active elems %d (%3.3f\%)",
1381 iter, nb_points, active_points, proc_nb_points,
1382 avtive_elems, proc_nb_active, avtive_full_elems,
1383 proc_nb_full_active, iter);
1384 }
1385 }
1386
1388 };
1389
1390 auto add_active_dofs_elem = [&](auto dm) {
1392 auto fe_pre_proc = boost::make_shared<FEMethod>();
1393 fe_pre_proc->preProcessHook = active_pre_lhs;
1394 auto fe_post_proc = boost::make_shared<FEMethod>();
1395 fe_post_proc->postProcessHook = active_post_lhs;
1397 ts_ctx_ptr->getPreProcessIJacobian().push_front(fe_pre_proc);
1398 ts_ctx_ptr->getPostProcessIJacobian().push_back(fe_post_proc);
1400 };
1401
1402 auto set_essential_bc = [&](auto dm, auto solver) {
1404
1405
1406 auto pre_proc_ptr = boost::make_shared<FEMethod>();
1407 auto post_proc_rhs_ptr = boost::make_shared<FEMethod>();
1408 auto post_proc_lhs_ptr = boost::make_shared<FEMethod>();
1410 ts_ctx_ptr->getPreProcessIFunction().push_front(pre_proc_ptr);
1411 ts_ctx_ptr->getPreProcessIJacobian().push_front(pre_proc_ptr);
1412 ts_ctx_ptr->getPostProcessIFunction().push_back(post_proc_rhs_ptr);
1413 ts_ctx_ptr->getPostProcessIJacobian().push_back(post_proc_lhs_ptr);
1414
1415
1416 auto disp_time_scale = boost::make_shared<TimeScale>();
1417
1418 auto get_bc_hook_rhs = [&]() {
1420 mField, pre_proc_ptr, {disp_time_scale},
false);
1421 };
1422 pre_proc_ptr->preProcessHook = get_bc_hook_rhs();
1423
1424 auto waak_post_proc_rhs_ptr = boost::weak_ptr<FEMethod>(
1425 post_proc_rhs_ptr);
1426 auto get_post_proc_hook_rhs = [this, waak_post_proc_rhs_ptr]() {
1429 mField, waak_post_proc_rhs_ptr.lock(),
nullptr, Sev::verbose)();
1431 mField, waak_post_proc_rhs_ptr.lock(), 1.)();
1433 };
1434 auto get_post_proc_hook_lhs = [&]() {
1436 mField, post_proc_lhs_ptr, 1.);
1437 };
1438
1439 post_proc_rhs_ptr->postProcessHook = get_post_proc_hook_rhs;
1440 post_proc_lhs_ptr->postProcessHook = get_post_proc_hook_lhs();
1441
1443 };
1444
1447 CHKERR TSSetIJacobian(solver,
B,
B, PETSC_NULLPTR, PETSC_NULLPTR);
1448 } else {
1449 CHKERR TSSetI2Jacobian(solver,
B,
B, PETSC_NULLPTR, PETSC_NULLPTR);
1450 }
1452 CHKERR TSSetSolution(solver,
D);
1453 } else {
1454 CHKERR TS2SetSolution(solver,
D, DD);
1455 }
1456 CHKERR set_section_monitor(solver);
1457 CHKERR set_time_monitor(dm, solver);
1458 CHKERR TSSetFromOptions(solver);
1459
1460 CHKERR add_active_dofs_elem(dm);
1461 boost::shared_ptr<SetUpSchur> schur_ptr;
1462 CHKERR set_schur_pc(solver, schur_ptr);
1463 CHKERR set_essential_bc(dm, solver);
1464
1468 BOOST_LOG_SCOPED_THREAD_ATTR("Timeline", attrs::timer());
1469 MOFEM_LOG(
"TIMER", Sev::verbose) <<
"TSSetUp";
1471 MOFEM_LOG(
"TIMER", Sev::verbose) <<
"TSSetUp <= done";
1472 MOFEM_LOG(
"TIMER", Sev::verbose) <<
"TSSolve";
1473 CHKERR TSSolve(solver, NULL);
1474 MOFEM_LOG(
"TIMER", Sev::verbose) <<
"TSSolve <= done";
1475
1479 "ts_manager_graph.dot");
1480 }
1481
1483}
#define CHK_MOAB_THROW(err, msg)
Check error code of MoAB function and throw MoFEM exception.
auto getDMTsCtx(DM dm)
Get TS context data structure used by DM.
MoFEMErrorCode MoFEMSNESMonitorFields(SNES snes, PetscInt its, PetscReal fgnorm, SnesCtx *ctx)
Sens monitor printing residual field by field.
auto getDMSubData(DM dm)
Get sub problem data structure.
boost::shared_ptr< BlockStructure > createBlockMatStructure(DM dm, SchurFEOpsFEandFields schur_fe_op_vec)
Create a Mat Diag Blocks object.
boost::shared_ptr< NestSchurData > createSchurNestedMatrixStruture(std::pair< SmartPetscObj< DM >, SmartPetscObj< DM > > dms, boost::shared_ptr< BlockStructure > block_mat_data_ptr, std::vector< std::string > fields_names, std::vector< boost::shared_ptr< Range > > field_ents, bool add_preconditioner_block)
Get the Schur Nest Mat Array object.
MoFEMErrorCode DMMoFEMSetNestSchurData(DM dm, boost::shared_ptr< NestSchurData >)
static auto getFTensor0FromVec(V &data)
Get tensor rank 0 (scalar) form data vector.
auto getDMSnesCtx(DM dm)
Get SNES context data structure used by DM.
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uYScatter
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uZScatter
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uXScatter
Section manager is used to create indexes and sections.
static MoFEMErrorCode writeTSGraphGraphviz(TsCtx *ts_ctx, std::string file_name)
TS graph to Graphviz file.
static std::array< int, 5 > activityData
static boost::shared_ptr< SetUpSchur > createSetUpSchur(MoFEM::Interface &m_field)
