833 {
835
839
841
842 auto set_section_monitor = [&](auto solver) {
844 SNES snes;
845 CHKERR TSGetSNES(solver, &snes);
846 CHKERR SNESMonitorSet(snes,
849 (void *)(snes_ctx_ptr.get()), nullptr);
851 };
852
853 auto create_post_process_elements = [&]() {
854 auto push_vol_ops = [this](auto &pip) {
856 pip, {
H1,
HDIV},
"GEOMETRY");
857
858 auto [common_plastic_ptr, common_hencky_ptr] =
859 PlasticOps::createCommonPlasticOps<SPACE_DIM, IT, DomainEleOp>(
860 mField,
"MAT_PLASTIC", pip,
"U",
"EP",
"TAU", 1., Sev::inform);
861
862 if (common_hencky_ptr) {
863 if (common_plastic_ptr->mGradPtr != common_hencky_ptr->matGradPtr)
865 }
866
867 return std::make_pair(common_plastic_ptr, common_hencky_ptr);
868 };
869
870 auto push_vol_post_proc_ops = [this](auto &pp_fe, auto &&p) {
872
873 auto &pip = pp_fe->getOpPtrVector();
874
875 auto [common_plastic_ptr, common_hencky_ptr] = p;
876
878
879 auto x_ptr = boost::make_shared<MatrixDouble>();
880 pip.push_back(
882 auto u_ptr = boost::make_shared<MatrixDouble>();
884
886
887 pip.push_back(
888
890
891 pp_fe->getPostProcMesh(), pp_fe->getMapGaussPts(),
892
893 {{"PLASTIC_SURFACE",
894 common_plastic_ptr->getPlasticSurfacePtr()},
895 {"PLASTIC_MULTIPLIER",
896 common_plastic_ptr->getPlasticTauPtr()}},
897
898 {{"U", u_ptr}, {"GEOMETRY", x_ptr}},
899
900 {{"GRAD", common_hencky_ptr->matGradPtr},
901 {"FIRST_PIOLA", common_hencky_ptr->getMatFirstPiolaStress()}},
902
903 {{"HENCKY_STRAIN", common_hencky_ptr->getMatLogC()},
904 {"PLASTIC_STRAIN", common_plastic_ptr->getPlasticStrainPtr()},
905 {"PLASTIC_FLOW", common_plastic_ptr->getPlasticFlowPtr()}}
906
907 )
908
909 );
910
911 } else {
912
913 pip.push_back(
914
916
917 pp_fe->getPostProcMesh(), pp_fe->getMapGaussPts(),
918
919 {{"PLASTIC_SURFACE",
920 common_plastic_ptr->getPlasticSurfacePtr()},
921 {"PLASTIC_MULTIPLIER",
922 common_plastic_ptr->getPlasticTauPtr()}},
923
924 {{"U", u_ptr}, {"GEOMETRY", x_ptr}},
925
926 {},
927
928 {{"STRAIN", common_plastic_ptr->mStrainPtr},
929 {"STRESS", common_plastic_ptr->mStressPtr},
930 {"PLASTIC_STRAIN", common_plastic_ptr->getPlasticStrainPtr()},
931 {"PLASTIC_FLOW", common_plastic_ptr->getPlasticFlowPtr()}}
932
933 )
934
935 );
936 }
937
939 };
940
941 PetscBool post_proc_vol;
942 PetscBool post_proc_skin;
943
945 post_proc_vol = PETSC_TRUE;
946 post_proc_skin = PETSC_FALSE;
947 } else {
948 post_proc_vol = PETSC_FALSE;
949 post_proc_skin = PETSC_TRUE;
950 }
952 PETSC_NULLPTR);
954 &post_proc_skin, PETSC_NULLPTR);
955
956 auto vol_post_proc = [this, push_vol_post_proc_ops, push_vol_ops,
957 post_proc_vol]() {
958 if (post_proc_vol == PETSC_FALSE)
959 return boost::shared_ptr<PostProcEle>();
960 auto pp_fe = boost::make_shared<PostProcEle>(
mField);
962 push_vol_post_proc_ops(pp_fe, push_vol_ops(pp_fe->getOpPtrVector())),
963 "push_vol_post_proc_ops");
964 return pp_fe;
965 };
966
967 auto skin_post_proc = [this, push_vol_post_proc_ops, push_vol_ops,
968 post_proc_skin]() {
969 if (post_proc_skin == PETSC_FALSE)
970 return boost::shared_ptr<SkinPostProcEle>();
971
973 auto pp_fe = boost::make_shared<SkinPostProcEle>(
mField);
976 pp_fe->getOpPtrVector().push_back(op_side);
978 pp_fe, push_vol_ops(op_side->getOpPtrVector())),
979 "push_vol_post_proc_ops");
980 return pp_fe;
981 };
982
983 return std::make_pair(vol_post_proc(), skin_post_proc());
984 };
985
986 auto scatter_create = [&](
auto D,
auto coeff) {
988 CHKERR is_manager->isCreateProblemFieldAndRank(
simple->getProblemName(),
989 ROW,
"U", coeff, coeff, is);
990 int loc_size;
991 CHKERR ISGetLocalSize(is, &loc_size);
994 VecScatter scatter;
995 CHKERR VecScatterCreate(
D, is,
v, PETSC_NULLPTR, &scatter);
998 };
999
1000 boost::shared_ptr<SetPtsData> field_eval_data;
1001 boost::shared_ptr<MatrixDouble> u_field_ptr;
1002
1003 std::array<double, 3> field_eval_coords{0.0, 0.0, 0.0};
1004 int coords_dim = 3;
1006 field_eval_coords.data(), &coords_dim,
1008
1009 boost::shared_ptr<std::map<std::string, boost::shared_ptr<VectorDouble>>>
1010 scalar_field_ptrs = boost::make_shared<
1011 std::map<std::string, boost::shared_ptr<VectorDouble>>>();
1012 boost::shared_ptr<std::map<std::string, boost::shared_ptr<MatrixDouble>>>
1013 vector_field_ptrs = boost::make_shared<
1014 std::map<std::string, boost::shared_ptr<MatrixDouble>>>();
1015 boost::shared_ptr<std::map<std::string, boost::shared_ptr<MatrixDouble>>>
1016 sym_tensor_field_ptrs = boost::make_shared<
1017 std::map<std::string, boost::shared_ptr<MatrixDouble>>>();
1018 boost::shared_ptr<std::map<std::string, boost::shared_ptr<MatrixDouble>>>
1019 tensor_field_ptrs = boost::make_shared<
1020 std::map<std::string, boost::shared_ptr<MatrixDouble>>>();
1021
1023 auto u_field_ptr = boost::make_shared<MatrixDouble>();
1024 field_eval_data =
1026
1028 field_eval_data,
simple->getDomainFEName());
1029
1030 field_eval_data->setEvalPoints(field_eval_coords.data(), 1);
1031 auto no_rule = [](
int,
int,
int) {
return -1; };
1032 auto field_eval_fe_ptr = field_eval_data->feMethodPtr;
1033 field_eval_fe_ptr->getRuleHook = no_rule;
1034
1036 field_eval_fe_ptr->getOpPtrVector(), {H1, HDIV}, "GEOMETRY");
1037
1038 auto [common_plastic_ptr, common_hencky_ptr] =
1039 PlasticOps::createCommonPlasticOps<SPACE_DIM, IT, DomainEleOp>(
1040 mField,
"MAT_PLASTIC", field_eval_fe_ptr->getOpPtrVector(),
"U",
1041 "EP", "TAU", 1., Sev::inform);
1042
1043 field_eval_fe_ptr->getOpPtrVector().push_back(
1045
1046 if ((common_plastic_ptr) && (common_hencky_ptr) && (scalar_field_ptrs)) {
1048 scalar_field_ptrs->insert(
1049 {"PLASTIC_SURFACE", common_plastic_ptr->getPlasticSurfacePtr()});
1050 scalar_field_ptrs->insert(
1051 {"PLASTIC_MULTIPLIER", common_plastic_ptr->getPlasticTauPtr()});
1052 vector_field_ptrs->insert({"U", u_field_ptr});
1053 sym_tensor_field_ptrs->insert(
1054 {"PLASTIC_STRAIN", common_plastic_ptr->getPlasticStrainPtr()});
1055 sym_tensor_field_ptrs->insert(
1056 {"PLASTIC_FLOW", common_plastic_ptr->getPlasticFlowPtr()});
1057 sym_tensor_field_ptrs->insert(
1058 {"HENCKY_STRAIN", common_hencky_ptr->getMatLogC()});
1059 tensor_field_ptrs->insert({"GRAD", common_hencky_ptr->matGradPtr});
1060 tensor_field_ptrs->insert(
1061 {"FIRST_PIOLA", common_hencky_ptr->getMatFirstPiolaStress()});
1062 } else {
1063 scalar_field_ptrs->insert(
1064 {"PLASTIC_SURFACE", common_plastic_ptr->getPlasticSurfacePtr()});
1065 scalar_field_ptrs->insert(
1066 {"PLASTIC_MULTIPLIER", common_plastic_ptr->getPlasticTauPtr()});
1067 vector_field_ptrs->insert({"U", u_field_ptr});
1068 sym_tensor_field_ptrs->insert(
1069 {"STRAIN", common_plastic_ptr->mStrainPtr});
1070 sym_tensor_field_ptrs->insert(
1071 {"STRESS", common_plastic_ptr->mStressPtr});
1072 sym_tensor_field_ptrs->insert(
1073 {"PLASTIC_STRAIN", common_plastic_ptr->getPlasticStrainPtr()});
1074 sym_tensor_field_ptrs->insert(
1075 {"PLASTIC_FLOW", common_plastic_ptr->getPlasticFlowPtr()});
1076 }
1077 }
1078 }
1079
1080 auto test_monitor_ptr = boost::make_shared<FEMethod>();
1081
1082 auto set_time_monitor = [&](auto dm, auto solver) {
1086 field_eval_coords, field_eval_data, scalar_field_ptrs,
1087 vector_field_ptrs, sym_tensor_field_ptrs, tensor_field_ptrs));
1088 boost::shared_ptr<ForcesAndSourcesCore> null;
1089
1090 test_monitor_ptr->postProcessHook = [&]() {
1092
1093 if (
atom_test && fabs(test_monitor_ptr->ts_t - 0.5) < 1e-12 &&
1094 test_monitor_ptr->ts_step == 25) {
1095
1096 if (scalar_field_ptrs->at("PLASTIC_MULTIPLIER")->size()) {
1097 auto t_tau =
1099 MOFEM_LOG(
"PlasticSync", Sev::inform) <<
"Eval point tau: " << t_tau;
1100
1101 if (
atom_test == 1 && fabs(t_tau - 0.688861) > 1e-5) {
1103 "atom test %d failed: wrong plastic multiplier value",
1105 }
1106 }
1107
1108 if (vector_field_ptrs->at("U")->size1()) {
1110 auto t_disp =
1111 getFTensor1FromMat<SPACE_DIM>(*vector_field_ptrs->at("U"));
1112 MOFEM_LOG(
"PlasticSync", Sev::inform) <<
"Eval point U: " << t_disp;
1113
1114 if (
atom_test == 1 && fabs(t_disp(0) - 0.25 / 2.) > 1e-5 ||
1115 fabs(t_disp(1) + 0.0526736) > 1e-5) {
1117 "atom test %d failed: wrong displacement value",
1119 }
1120 }
1121
1122 if (sym_tensor_field_ptrs->at("PLASTIC_STRAIN")->size1()) {
1123 auto t_plastic_strain = getFTensor2SymmetricFromMat<SPACE_DIM>(
1124 *sym_tensor_field_ptrs->at("PLASTIC_STRAIN"));
1126 << "Eval point EP: " << t_plastic_strain;
1127
1129 fabs(t_plastic_strain(0, 0) - 0.221943) > 1e-5 ||
1130 fabs(t_plastic_strain(0, 1)) > 1e-5 ||
1131 fabs(t_plastic_strain(1, 1) + 0.110971) > 1e-5) {
1133 "atom test %d failed: wrong plastic strain value",
1135 }
1136 }
1137
1138 if (tensor_field_ptrs->at("FIRST_PIOLA")->size1()) {
1139 auto t_piola_stress = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(
1140 *tensor_field_ptrs->at("FIRST_PIOLA"));
1142 << "Eval point Piola stress: " << t_piola_stress;
1143
1144 if (
atom_test == 1 && fabs((t_piola_stress(0, 0) - 198.775) /
1145 t_piola_stress(0, 0)) > 1e-5 ||
1146 fabs(t_piola_stress(0, 1)) + fabs(t_piola_stress(1, 0)) +
1147 fabs(t_piola_stress(1, 1)) >
1148 1e-5) {
1150 "atom test %d failed: wrong Piola stress value",
1152 }
1153 }
1154 }
1155
1158 };
1159
1161 monitor_ptr, null, test_monitor_ptr);
1162
1164 };
1165
1166 auto set_schur_pc = [&](auto solver,
1167 boost::shared_ptr<SetUpSchur> &schur_ptr) {
1169
1170 auto name_prb =
simple->getProblemName();
1171
1172
1181 for (
auto f : {
"U"}) {
1184 }
1186
1188 };
1189
1198#ifdef ADD_CONTACT
1199 for (
auto f : {
"SIGMA",
"EP",
"TAU"}) {
1202 }
1203#else
1204 for (
auto f : {
"EP",
"TAU"}) {
1207 }
1208#endif
1211 };
1212
1213
1214 if constexpr (
AT == AssemblyType::BLOCK_SCHUR) {
1215
1220
1221#ifdef ADD_CONTACT
1222
1223 auto get_nested_mat_data = [&](auto schur_dm, auto block_dm) {
1226
1227 {
1228
1229 {simple->getDomainFEName(),
1230
1231 {{"U", "U"},
1232 {"SIGMA", "SIGMA"},
1233 {"U", "SIGMA"},
1234 {"SIGMA", "U"},
1235 {"EP", "EP"},
1236 {"TAU", "TAU"},
1237 {"U", "EP"},
1238 {"EP", "U"},
1239 {"EP", "TAU"},
1240 {"TAU", "EP"},
1241 {"TAU", "U"}
1242
1243 }},
1244
1245 {simple->getBoundaryFEName(),
1246
1247 {{"SIGMA", "SIGMA"}, {"U", "SIGMA"}, {"SIGMA", "U"}
1248
1249 }}
1250
1251 }
1252
1253 );
1254
1256
1257 {dm_schur, dm_block}, block_mat_data,
1258
1259 {"SIGMA", "EP", "TAU"}, {nullptr, nullptr, nullptr}, true
1260
1261 );
1262 };
1263
1264#else
1265
1266 auto get_nested_mat_data = [&](auto schur_dm, auto block_dm) {
1267 auto block_mat_data =
1269
1270 {{simple->getDomainFEName(),
1271
1272 {{"U", "U"},
1273 {"EP", "EP"},
1274 {"TAU", "TAU"},
1275 {"U", "EP"},
1276 {"EP", "U"},
1277 {"EP", "TAU"},
1278 {"TAU", "U"},
1279 {"TAU", "EP"}
1280
1281 }}}
1282
1283 );
1284
1286
1287 {dm_schur, dm_block}, block_mat_data,
1288
1289 {"EP", "TAU"}, {nullptr, nullptr}, false
1290
1291 );
1292 };
1293
1294#endif
1295
1296 auto nested_mat_data = get_nested_mat_data(dm_schur, dm_block);
1298
1299 auto block_is =
getDMSubData(dm_block)->getSmartRowIs();
1300 auto ao_schur =
getDMSubData(dm_schur)->getSmartRowMap();
1301
1302
1303
1304 schur_ptr =
1306 CHKERR schur_ptr->setUp(solver);
1307 }
1308
1310 };
1311
1312 auto dm =
simple->getDM();
1315 CHKERR VecSetDM(
D, PETSC_NULLPTR);
1316 CHKERR VecSetDM(DD, PETSC_NULLPTR);
1321
1322 auto create_solver = [pip_mng]() {
1324 return pip_mng->createTSIM();
1325 else
1326 return pip_mng->createTSIM2();
1327 };
1328
1329 auto solver = create_solver();
1330
1331 auto active_pre_lhs = []() {
1336 };
1337
1338 auto active_post_lhs = [&]() {
1340 auto get_iter = [&]() {
1341 SNES snes;
1343 int iter;
1345 "Can not get iter");
1346 return iter;
1347 };
1348
1349 auto iter = get_iter();
1350 if (iter >= 0) {
1351
1352 std::array<int, 5> activity_data;
1353 std::fill(activity_data.begin(), activity_data.end(), 0);
1355 activity_data.data(), activity_data.size(), MPI_INT,
1357
1358 int &active_points = activity_data[0];
1359 int &avtive_full_elems = activity_data[1];
1360 int &avtive_elems = activity_data[2];
1361 int &nb_points = activity_data[3];
1362 int &nb_elements = activity_data[4];
1363
1364 if (nb_points) {
1365
1366 double proc_nb_points =
1367 100 * static_cast<double>(active_points) / nb_points;
1368 double proc_nb_active =
1369 100 * static_cast<double>(avtive_elems) / nb_elements;
1370 double proc_nb_full_active = 100;
1371 if (avtive_elems)
1372 proc_nb_full_active =
1373 100 * static_cast<double>(avtive_full_elems) / avtive_elems;
1374
1376 "Iter %d nb pts %d nb active pts %d (%3.3f\%) nb active "
1377 "elements %d "
1378 "(%3.3f\%) nb full active elems %d (%3.3f\%)",
1379 iter, nb_points, active_points, proc_nb_points,
1380 avtive_elems, proc_nb_active, avtive_full_elems,
1381 proc_nb_full_active, iter);
1382 }
1383 }
1384
1386 };
1387
1388 auto add_active_dofs_elem = [&](auto dm) {
1390 auto fe_pre_proc = boost::make_shared<FEMethod>();
1391 fe_pre_proc->preProcessHook = active_pre_lhs;
1392 auto fe_post_proc = boost::make_shared<FEMethod>();
1393 fe_post_proc->postProcessHook = active_post_lhs;
1395 ts_ctx_ptr->getPreProcessIJacobian().push_front(fe_pre_proc);
1396 ts_ctx_ptr->getPostProcessIJacobian().push_back(fe_post_proc);
1398 };
1399
1400 auto set_essential_bc = [&](auto dm, auto solver) {
1402
1403
1404 auto pre_proc_ptr = boost::make_shared<FEMethod>();
1405 auto post_proc_rhs_ptr = boost::make_shared<FEMethod>();
1406 auto post_proc_lhs_ptr = boost::make_shared<FEMethod>();
1408 ts_ctx_ptr->getPreProcessIFunction().push_front(pre_proc_ptr);
1409 ts_ctx_ptr->getPreProcessIJacobian().push_front(pre_proc_ptr);
1410 ts_ctx_ptr->getPostProcessIFunction().push_back(post_proc_rhs_ptr);
1411 ts_ctx_ptr->getPostProcessIJacobian().push_back(post_proc_lhs_ptr);
1412
1413
1414 auto disp_time_scale = boost::make_shared<TimeScale>();
1415
1416 auto get_bc_hook_rhs = [&]() {
1418 mField, pre_proc_ptr, {disp_time_scale},
false);
1419 };
1420 pre_proc_ptr->preProcessHook = get_bc_hook_rhs();
1421
1422 auto waak_post_proc_rhs_ptr = boost::weak_ptr<FEMethod>(
1423 post_proc_rhs_ptr);
1424 auto get_post_proc_hook_rhs = [this, waak_post_proc_rhs_ptr]() {
1427 mField, waak_post_proc_rhs_ptr.lock(),
nullptr, Sev::verbose)();
1429 mField, waak_post_proc_rhs_ptr.lock(), 1.)();
1431 };
1432 auto get_post_proc_hook_lhs = [&]() {
1434 mField, post_proc_lhs_ptr, 1.);
1435 };
1436
1437 post_proc_rhs_ptr->postProcessHook = get_post_proc_hook_rhs;
1438 post_proc_lhs_ptr->postProcessHook = get_post_proc_hook_lhs();
1439
1441 };
1442
1445 CHKERR TSSetIJacobian(solver,
B,
B, PETSC_NULLPTR, PETSC_NULLPTR);
1446 } else {
1447 CHKERR TSSetI2Jacobian(solver,
B,
B, PETSC_NULLPTR, PETSC_NULLPTR);
1448 }
1450 CHKERR TSSetSolution(solver,
D);
1451 } else {
1452 CHKERR TS2SetSolution(solver,
D, DD);
1453 }
1454 CHKERR set_section_monitor(solver);
1455 CHKERR set_time_monitor(dm, solver);
1456 CHKERR TSSetFromOptions(solver);
1457
1458 CHKERR add_active_dofs_elem(dm);
1459 boost::shared_ptr<SetUpSchur> schur_ptr;
1460 CHKERR set_schur_pc(solver, schur_ptr);
1461 CHKERR set_essential_bc(dm, solver);
1462
1466 BOOST_LOG_SCOPED_THREAD_ATTR("Timeline", attrs::timer());
1467 MOFEM_LOG(
"TIMER", Sev::verbose) <<
"TSSetUp";
1469 MOFEM_LOG(
"TIMER", Sev::verbose) <<
"TSSetUp <= done";
1470 MOFEM_LOG(
"TIMER", Sev::verbose) <<
"TSSolve";
1471 CHKERR TSSolve(solver, NULL);
1472 MOFEM_LOG(
"TIMER", Sev::verbose) <<
"TSSolve <= done";
1473
1477 "ts_manager_graph.dot");
1478 }
1479
1481}
#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)
