834 {
836
840
842
843 auto set_section_monitor = [&](auto solver) {
845 SNES snes;
846 CHKERR TSGetSNES(solver, &snes);
847 CHKERR SNESMonitorSet(snes,
850 (void *)(snes_ctx_ptr.get()), nullptr);
852 };
853
854 auto create_post_process_elements = [&]() {
855 auto push_vol_ops = [this](auto &pip) {
857 pip, {
H1,
HDIV},
"GEOMETRY");
858
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);
862
863 if (common_hencky_ptr) {
864 if (common_plastic_ptr->mGradPtr != common_hencky_ptr->matGradPtr)
866 }
867
868 return std::make_pair(common_plastic_ptr, common_hencky_ptr);
869 };
870
871 auto push_vol_post_proc_ops = [this](auto &pp_fe, auto &&p) {
873
874 auto &pip = pp_fe->getOpPtrVector();
875
876 auto [common_plastic_ptr, common_hencky_ptr] = p;
877
879
880 auto x_ptr = boost::make_shared<MatrixDouble>();
881 pip.push_back(
883 auto u_ptr = boost::make_shared<MatrixDouble>();
885
887
888 pip.push_back(
889
891
892 pp_fe->getPostProcMesh(), pp_fe->getMapGaussPts(),
893
894 {{"PLASTIC_SURFACE",
895 common_plastic_ptr->getPlasticSurfacePtr()},
896 {"PLASTIC_MULTIPLIER",
897 common_plastic_ptr->getPlasticTauPtr()}},
898
899 {{"U", u_ptr}, {"GEOMETRY", x_ptr}},
900
901 {{"GRAD", common_hencky_ptr->matGradPtr},
902 {"FIRST_PIOLA", common_hencky_ptr->getMatFirstPiolaStress()}},
903
904 {{"HENCKY_STRAIN", common_hencky_ptr->getMatLogC()},
905 {"PLASTIC_STRAIN", common_plastic_ptr->getPlasticStrainPtr()},
906 {"PLASTIC_FLOW", common_plastic_ptr->getPlasticFlowPtr()}}
907
908 )
909
910 );
911
912 } else {
913
914 pip.push_back(
915
917
918 pp_fe->getPostProcMesh(), pp_fe->getMapGaussPts(),
919
920 {{"PLASTIC_SURFACE",
921 common_plastic_ptr->getPlasticSurfacePtr()},
922 {"PLASTIC_MULTIPLIER",
923 common_plastic_ptr->getPlasticTauPtr()}},
924
925 {{"U", u_ptr}, {"GEOMETRY", x_ptr}},
926
927 {},
928
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()}}
933
934 )
935
936 );
937 }
938
940 };
941
942 PetscBool post_proc_vol;
943 PetscBool post_proc_skin;
944
946 post_proc_vol = PETSC_TRUE;
947 post_proc_skin = PETSC_FALSE;
948 } else {
949 post_proc_vol = PETSC_FALSE;
950 post_proc_skin = PETSC_TRUE;
951 }
953 PETSC_NULL);
955 &post_proc_skin, PETSC_NULL);
956
957 auto vol_post_proc = [this, push_vol_post_proc_ops, push_vol_ops,
958 post_proc_vol]() {
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");
965 return pp_fe;
966 };
967
968 auto skin_post_proc = [this, push_vol_post_proc_ops, push_vol_ops,
969 post_proc_skin]() {
970 if (post_proc_skin == PETSC_FALSE)
971 return boost::shared_ptr<SkinPostProcEle>();
972
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");
981 return pp_fe;
982 };
983
984 return std::make_pair(vol_post_proc(), skin_post_proc());
985 };
986
987 auto scatter_create = [&](
auto D,
auto coeff) {
989 CHKERR is_manager->isCreateProblemFieldAndRank(
simple->getProblemName(),
990 ROW,
"U", coeff, coeff, is);
991 int loc_size;
992 CHKERR ISGetLocalSize(is, &loc_size);
995 VecScatter scatter;
996 CHKERR VecScatterCreate(
D, is,
v, PETSC_NULL, &scatter);
999 };
1000
1001 boost::shared_ptr<SetPtsData> field_eval_data;
1002 boost::shared_ptr<MatrixDouble> u_field_ptr;
1003
1004 std::array<double, SPACE_DIM> field_eval_coords;
1007 field_eval_coords.data(), &coords_dim,
1009
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>>>();
1022
1024 auto u_field_ptr = boost::make_shared<MatrixDouble>();
1025 field_eval_data =
1027
1029 field_eval_data,
simple->getDomainFEName());
1030
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;
1035
1037 field_eval_fe_ptr->getOpPtrVector(), {H1, HDIV}, "GEOMETRY");
1038
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);
1043
1044 field_eval_fe_ptr->getOpPtrVector().push_back(
1046
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()});
1063 } else {
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()});
1077 }
1078 }
1079 }
1080
1081 auto test_monitor_ptr = boost::make_shared<FEMethod>();
1082
1083 auto set_time_monitor = [&](auto dm, auto solver) {
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;
1090
1091 test_monitor_ptr->postProcessHook = [&]() {
1093
1094 if (
atom_test && fabs(test_monitor_ptr->ts_t - 0.5) < 1e-12 &&
1095 test_monitor_ptr->ts_step == 25) {
1096
1097 if (scalar_field_ptrs->at("PLASTIC_MULTIPLIER")->size()) {
1098 auto t_tau =
1100 MOFEM_LOG(
"PlasticSync", Sev::inform) <<
"Eval point tau: " << t_tau;
1101
1102 if (
atom_test == 1 && fabs(t_tau - 0.688861) > 1e-5) {
1104 "atom test %d failed: wrong plastic multiplier value",
1106 }
1107 }
1108
1109 if (vector_field_ptrs->at("U")->size1()) {
1111 auto t_disp =
1112 getFTensor1FromMat<SPACE_DIM>(*vector_field_ptrs->at("U"));
1113 MOFEM_LOG(
"PlasticSync", Sev::inform) <<
"Eval point U: " << t_disp;
1114
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",
1120 }
1121 }
1122
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;
1128
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",
1136 }
1137 }
1138
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;
1144
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)) >
1149 1e-5) {
1151 "atom test %d failed: wrong Piola stress value",
1153 }
1154 }
1155 }
1156
1159 };
1160
1162 monitor_ptr, null, test_monitor_ptr);
1163
1165 };
1166
1167 auto set_schur_pc = [&](auto solver,
1168 boost::shared_ptr<SetUpSchur> &schur_ptr) {
1170
1171 auto name_prb =
simple->getProblemName();
1172
1173
1182 for (
auto f : {
"U"}) {
1185 }
1187
1189 };
1190
1199#ifdef ADD_CONTACT
1200 for (
auto f : {
"SIGMA",
"EP",
"TAU"}) {
1203 }
1204#else
1205 for (
auto f : {
"EP",
"TAU"}) {
1208 }
1209#endif
1212 };
1213
1214
1215 if constexpr (
AT == AssemblyType::BLOCK_SCHUR) {
1216
1221
1222#ifdef ADD_CONTACT
1223
1224 auto get_nested_mat_data = [&](auto schur_dm, auto block_dm) {
1227
1228 {
1229
1230 {simple->getDomainFEName(),
1231
1232 {{"U", "U"},
1233 {"SIGMA", "SIGMA"},
1234 {"U", "SIGMA"},
1235 {"SIGMA", "U"},
1236 {"EP", "EP"},
1237 {"TAU", "TAU"},
1238 {"U", "EP"},
1239 {"EP", "U"},
1240 {"EP", "TAU"},
1241 {"TAU", "EP"},
1242 {"TAU", "U"}
1243
1244 }},
1245
1246 {simple->getBoundaryFEName(),
1247
1248 {{"SIGMA", "SIGMA"}, {"U", "SIGMA"}, {"SIGMA", "U"}
1249
1250 }}
1251
1252 }
1253
1254 );
1255
1257
1258 {dm_schur, dm_block}, block_mat_data,
1259
1260 {"SIGMA", "EP", "TAU"}, {nullptr, nullptr, nullptr}, true
1261
1262 );
1263 };
1264
1265#else
1266
1267 auto get_nested_mat_data = [&](auto schur_dm, auto block_dm) {
1268 auto block_mat_data =
1270
1271 {{simple->getDomainFEName(),
1272
1273 {{"U", "U"},
1274 {"EP", "EP"},
1275 {"TAU", "TAU"},
1276 {"U", "EP"},
1277 {"EP", "U"},
1278 {"EP", "TAU"},
1279 {"TAU", "U"},
1280 {"TAU", "EP"}
1281
1282 }}}
1283
1284 );
1285
1287
1288 {dm_schur, dm_block}, block_mat_data,
1289
1290 {"EP", "TAU"}, {nullptr, nullptr}, false
1291
1292 );
1293 };
1294
1295#endif
1296
1297 auto nested_mat_data = get_nested_mat_data(dm_schur, dm_block);
1299
1300 auto block_is =
getDMSubData(dm_block)->getSmartRowIs();
1301 auto ao_schur =
getDMSubData(dm_schur)->getSmartRowMap();
1302
1303
1304
1305 schur_ptr =
1307 CHKERR schur_ptr->setUp(solver);
1308 }
1309
1311 };
1312
1313 auto dm =
simple->getDM();
1320
1321 auto create_solver = [pip_mng]() {
1323 return pip_mng->createTSIM();
1324 else
1325 return pip_mng->createTSIM2();
1326 };
1327
1328 auto solver = create_solver();
1329
1330 auto active_pre_lhs = []() {
1335 };
1336
1337 auto active_post_lhs = [&]() {
1339 auto get_iter = [&]() {
1340 SNES snes;
1342 int iter;
1344 "Can not get iter");
1345 return iter;
1346 };
1347
1348 auto iter = get_iter();
1349 if (iter >= 0) {
1350
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,
1356
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];
1362
1363 if (nb_points) {
1364
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;
1370 if (avtive_elems)
1371 proc_nb_full_active =
1372 100 * static_cast<double>(avtive_full_elems) / avtive_elems;
1373
1375 "Iter %d nb pts %d nb active pts %d (%3.3f\%) nb active "
1376 "elements %d "
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);
1381 }
1382 }
1383
1385 };
1386
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);
1397 };
1398
1399 auto set_essential_bc = [&](auto dm, auto solver) {
1401
1402
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>();
1406
1407
1408 auto disp_time_scale = boost::make_shared<TimeScale>();
1409
1410 auto get_bc_hook_rhs = [this, pre_proc_ptr, disp_time_scale]() {
1412 {disp_time_scale}, false);
1413 return hook;
1414 };
1415 pre_proc_ptr->preProcessHook = get_bc_hook_rhs();
1416
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.)();
1424 };
1425 auto get_post_proc_hook_lhs = [this, post_proc_lhs_ptr]() {
1427 mField, post_proc_lhs_ptr, 1.);
1428 };
1429 post_proc_rhs_ptr->postProcessHook = get_post_proc_hook_rhs;
1430
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);
1437
1439 };
1440
1442 CHKERR TSSetSolution(solver,
D);
1443 } else {
1444 CHKERR TS2SetSolution(solver,
D, DD);
1445 }
1446
1447 CHKERR set_section_monitor(solver);
1448 CHKERR set_time_monitor(dm, solver);
1449 CHKERR TSSetFromOptions(solver);
1451
1452 CHKERR add_active_dofs_elem(dm);
1453 boost::shared_ptr<SetUpSchur> schur_ptr;
1454 CHKERR set_schur_pc(solver, schur_ptr);
1455 CHKERR set_essential_bc(dm, solver);
1456
1460 BOOST_LOG_SCOPED_THREAD_ATTR("Timeline", attrs::timer());
1461 MOFEM_LOG(
"TIMER", Sev::verbose) <<
"TSSetUp";
1463 MOFEM_LOG(
"TIMER", Sev::verbose) <<
"TSSetUp <= done";
1464 MOFEM_LOG(
"TIMER", Sev::verbose) <<
"TSSolve";
1465 CHKERR TSSolve(solver, NULL);
1466 MOFEM_LOG(
"TIMER", Sev::verbose) <<
"TSSolve <= done";
1467
1469}
1470
1471
1472
#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.
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
MoFEMErrorCode DMMoFEMSetNestSchurData(DM dm, boost::shared_ptr< NestSchurData >)
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 std::array< int, 5 > activityData
static boost::shared_ptr< SetUpSchur > createSetUpSchur(MoFEM::Interface &m_field)
