832 {
834
838
840
841 auto set_section_monitor = [&](auto solver) {
843 SNES snes;
844 CHKERR TSGetSNES(solver, &snes);
845 CHKERR SNESMonitorSet(snes,
848 (void *)(snes_ctx_ptr.get()), nullptr);
850 };
851
852 auto create_post_process_elements = [&]() {
853 auto push_vol_ops = [this](auto &pip) {
855 pip, {
H1,
HDIV},
"GEOMETRY");
856
857 auto [common_plastic_ptr, common_hencky_ptr] =
858 PlasticOps::createCommonPlasticOps<SPACE_DIM, IT, DomainEleOp>(
859 mField,
"MAT_PLASTIC", pip,
"U",
"EP",
"TAU", 1., Sev::inform);
860
861 if (common_hencky_ptr) {
862 if (common_plastic_ptr->mGradPtr != common_hencky_ptr->matGradPtr)
864 }
865
866 return std::make_pair(common_plastic_ptr, common_hencky_ptr);
867 };
868
869 auto push_vol_post_proc_ops = [this](auto &pp_fe, auto &&p) {
871
872 auto &pip = pp_fe->getOpPtrVector();
873
874 auto [common_plastic_ptr, common_hencky_ptr] = p;
875
877
878 auto x_ptr = boost::make_shared<MatrixDouble>();
879 pip.push_back(
881 auto u_ptr = boost::make_shared<MatrixDouble>();
883
885
886 pip.push_back(
887
889
890 pp_fe->getPostProcMesh(), pp_fe->getMapGaussPts(),
891
892 {{"PLASTIC_SURFACE",
893 common_plastic_ptr->getPlasticSurfacePtr()},
894 {"PLASTIC_MULTIPLIER",
895 common_plastic_ptr->getPlasticTauPtr()}},
896
897 {{"U", u_ptr}, {"GEOMETRY", x_ptr}},
898
899 {{"GRAD", common_hencky_ptr->matGradPtr},
900 {"FIRST_PIOLA", common_hencky_ptr->getMatFirstPiolaStress()}},
901
902 {{"HENCKY_STRAIN", common_hencky_ptr->getMatLogC()},
903 {"PLASTIC_STRAIN", common_plastic_ptr->getPlasticStrainPtr()},
904 {"PLASTIC_FLOW", common_plastic_ptr->getPlasticFlowPtr()}}
905
906 )
907
908 );
909
910 } else {
911
912 pip.push_back(
913
915
916 pp_fe->getPostProcMesh(), pp_fe->getMapGaussPts(),
917
918 {{"PLASTIC_SURFACE",
919 common_plastic_ptr->getPlasticSurfacePtr()},
920 {"PLASTIC_MULTIPLIER",
921 common_plastic_ptr->getPlasticTauPtr()}},
922
923 {{"U", u_ptr}, {"GEOMETRY", x_ptr}},
924
925 {},
926
927 {{"STRAIN", common_plastic_ptr->mStrainPtr},
928 {"STRESS", common_plastic_ptr->mStressPtr},
929 {"PLASTIC_STRAIN", common_plastic_ptr->getPlasticStrainPtr()},
930 {"PLASTIC_FLOW", common_plastic_ptr->getPlasticFlowPtr()}}
931
932 )
933
934 );
935 }
936
938 };
939
940 PetscBool post_proc_vol;
941 PetscBool post_proc_skin;
942
944 post_proc_vol = PETSC_TRUE;
945 post_proc_skin = PETSC_FALSE;
946 } else {
947 post_proc_vol = PETSC_FALSE;
948 post_proc_skin = PETSC_TRUE;
949 }
951 PETSC_NULLPTR);
953 &post_proc_skin, PETSC_NULLPTR);
954
955 auto vol_post_proc = [this, push_vol_post_proc_ops, push_vol_ops,
956 post_proc_vol]() {
957 if (post_proc_vol == PETSC_FALSE)
958 return boost::shared_ptr<PostProcEle>();
959 auto pp_fe = boost::make_shared<PostProcEle>(
mField);
961 push_vol_post_proc_ops(pp_fe, push_vol_ops(pp_fe->getOpPtrVector())),
962 "push_vol_post_proc_ops");
963 return pp_fe;
964 };
965
966 auto skin_post_proc = [this, push_vol_post_proc_ops, push_vol_ops,
967 post_proc_skin]() {
968 if (post_proc_skin == PETSC_FALSE)
969 return boost::shared_ptr<SkinPostProcEle>();
970
972 auto pp_fe = boost::make_shared<SkinPostProcEle>(
mField);
975 pp_fe->getOpPtrVector().push_back(op_side);
977 pp_fe, push_vol_ops(op_side->getOpPtrVector())),
978 "push_vol_post_proc_ops");
979 return pp_fe;
980 };
981
982 return std::make_pair(vol_post_proc(), skin_post_proc());
983 };
984
985 auto scatter_create = [&](
auto D,
auto coeff) {
988 ROW,
"U", coeff, coeff, is);
989 int loc_size;
990 CHKERR ISGetLocalSize(is, &loc_size);
993 VecScatter scatter;
994 CHKERR VecScatterCreate(
D, is,
v, PETSC_NULLPTR, &scatter);
997 };
998
999 boost::shared_ptr<SetPtsData> field_eval_data;
1000 boost::shared_ptr<MatrixDouble> u_field_ptr;
1001
1002 std::array<double, 3> field_eval_coords{0.0, 0.0, 0.0};
1003 int coords_dim = 3;
1005 field_eval_coords.data(), &coords_dim,
1007
1008 boost::shared_ptr<std::map<std::string, boost::shared_ptr<VectorDouble>>>
1009 scalar_field_ptrs = boost::make_shared<
1010 std::map<std::string, boost::shared_ptr<VectorDouble>>>();
1011 boost::shared_ptr<std::map<std::string, boost::shared_ptr<MatrixDouble>>>
1012 vector_field_ptrs = boost::make_shared<
1013 std::map<std::string, boost::shared_ptr<MatrixDouble>>>();
1014 boost::shared_ptr<std::map<std::string, boost::shared_ptr<MatrixDouble>>>
1015 sym_tensor_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 tensor_field_ptrs = boost::make_shared<
1019 std::map<std::string, boost::shared_ptr<MatrixDouble>>>();
1020
1022 auto u_field_ptr = boost::make_shared<MatrixDouble>();
1023 field_eval_data =
1025
1028
1029 field_eval_data->setEvalPoints(field_eval_coords.data(), 1);
1030 auto no_rule = [](
int,
int,
int) {
return -1; };
1031 auto field_eval_fe_ptr = field_eval_data->feMethodPtr;
1032 field_eval_fe_ptr->getRuleHook = no_rule;
1033
1035 field_eval_fe_ptr->getOpPtrVector(), {H1, HDIV}, "GEOMETRY");
1036
1037 auto [common_plastic_ptr, common_hencky_ptr] =
1038 PlasticOps::createCommonPlasticOps<SPACE_DIM, IT, DomainEleOp>(
1039 mField,
"MAT_PLASTIC", field_eval_fe_ptr->getOpPtrVector(),
"U",
1040 "EP", "TAU", 1., Sev::inform);
1041
1042 field_eval_fe_ptr->getOpPtrVector().push_back(
1044
1045 if ((common_plastic_ptr) && (common_hencky_ptr) && (scalar_field_ptrs)) {
1047 scalar_field_ptrs->insert(
1048 {"PLASTIC_SURFACE", common_plastic_ptr->getPlasticSurfacePtr()});
1049 scalar_field_ptrs->insert(
1050 {"PLASTIC_MULTIPLIER", common_plastic_ptr->getPlasticTauPtr()});
1051 vector_field_ptrs->insert({"U", u_field_ptr});
1052 sym_tensor_field_ptrs->insert(
1053 {"PLASTIC_STRAIN", common_plastic_ptr->getPlasticStrainPtr()});
1054 sym_tensor_field_ptrs->insert(
1055 {"PLASTIC_FLOW", common_plastic_ptr->getPlasticFlowPtr()});
1056 sym_tensor_field_ptrs->insert(
1057 {"HENCKY_STRAIN", common_hencky_ptr->getMatLogC()});
1058 tensor_field_ptrs->insert({"GRAD", common_hencky_ptr->matGradPtr});
1059 tensor_field_ptrs->insert(
1060 {"FIRST_PIOLA", common_hencky_ptr->getMatFirstPiolaStress()});
1061 } else {
1062 scalar_field_ptrs->insert(
1063 {"PLASTIC_SURFACE", common_plastic_ptr->getPlasticSurfacePtr()});
1064 scalar_field_ptrs->insert(
1065 {"PLASTIC_MULTIPLIER", common_plastic_ptr->getPlasticTauPtr()});
1066 vector_field_ptrs->insert({"U", u_field_ptr});
1067 sym_tensor_field_ptrs->insert(
1068 {"STRAIN", common_plastic_ptr->mStrainPtr});
1069 sym_tensor_field_ptrs->insert(
1070 {"STRESS", common_plastic_ptr->mStressPtr});
1071 sym_tensor_field_ptrs->insert(
1072 {"PLASTIC_STRAIN", common_plastic_ptr->getPlasticStrainPtr()});
1073 sym_tensor_field_ptrs->insert(
1074 {"PLASTIC_FLOW", common_plastic_ptr->getPlasticFlowPtr()});
1075 }
1076 }
1077 }
1078
1079 auto test_monitor_ptr = boost::make_shared<FEMethod>();
1080
1081 auto set_time_monitor = [&](auto dm, auto solver) {
1085 field_eval_coords, field_eval_data, scalar_field_ptrs,
1086 vector_field_ptrs, sym_tensor_field_ptrs, tensor_field_ptrs));
1087 boost::shared_ptr<ForcesAndSourcesCore> null;
1088
1089 test_monitor_ptr->postProcessHook = [&]() {
1091
1092 if (
atom_test && fabs(test_monitor_ptr->ts_t - 0.5) < 1e-12 &&
1093 test_monitor_ptr->ts_step == 25) {
1094
1095 if (scalar_field_ptrs->at("PLASTIC_MULTIPLIER")->size()) {
1096 auto t_tau =
1098 MOFEM_LOG(
"PlasticSync", Sev::inform) <<
"Eval point tau: " << t_tau;
1099
1100 if (
atom_test == 1 && fabs(t_tau - 0.688861) > 1e-5) {
1102 "atom test %d failed: wrong plastic multiplier value",
1104 }
1105 }
1106
1107 if (vector_field_ptrs->at("U")->size1()) {
1109 auto t_disp =
1110 getFTensor1FromMat<SPACE_DIM>(*vector_field_ptrs->at("U"));
1111 MOFEM_LOG(
"PlasticSync", Sev::inform) <<
"Eval point U: " << t_disp;
1112
1113 if (
atom_test == 1 && fabs(t_disp(0) - 0.25 / 2.) > 1e-5 ||
1114 fabs(t_disp(1) + 0.0526736) > 1e-5) {
1116 "atom test %d failed: wrong displacement value",
1118 }
1119 }
1120
1121 if (sym_tensor_field_ptrs->at("PLASTIC_STRAIN")->size1()) {
1122 auto t_plastic_strain = getFTensor2SymmetricFromMat<SPACE_DIM>(
1123 *sym_tensor_field_ptrs->at("PLASTIC_STRAIN"));
1125 << "Eval point EP: " << t_plastic_strain;
1126
1128 fabs(t_plastic_strain(0, 0) - 0.221943) > 1e-5 ||
1129 fabs(t_plastic_strain(0, 1)) > 1e-5 ||
1130 fabs(t_plastic_strain(1, 1) + 0.110971) > 1e-5) {
1132 "atom test %d failed: wrong plastic strain value",
1134 }
1135 }
1136
1137 if (tensor_field_ptrs->at("FIRST_PIOLA")->size1()) {
1138 auto t_piola_stress = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(
1139 *tensor_field_ptrs->at("FIRST_PIOLA"));
1141 << "Eval point Piola stress: " << t_piola_stress;
1142
1143 if (
atom_test == 1 && fabs((t_piola_stress(0, 0) - 198.775) /
1144 t_piola_stress(0, 0)) > 1e-5 ||
1145 fabs(t_piola_stress(0, 1)) + fabs(t_piola_stress(1, 0)) +
1146 fabs(t_piola_stress(1, 1)) >
1147 1e-5) {
1149 "atom test %d failed: wrong Piola stress value",
1151 }
1152 }
1153 }
1154
1157 };
1158
1160 monitor_ptr, null, test_monitor_ptr);
1161
1163 };
1164
1165 auto set_schur_pc = [&](auto solver,
1166 boost::shared_ptr<SetUpSchur> &schur_ptr) {
1168
1170
1171
1180 for (
auto f : {
"U"}) {
1183 }
1185
1187 };
1188
1197#ifdef ADD_CONTACT
1198 for (
auto f : {
"SIGMA",
"EP",
"TAU"}) {
1201 }
1202#else
1203 for (
auto f : {
"EP",
"TAU"}) {
1206 }
1207#endif
1210 };
1211
1212
1213 if constexpr (
AT == AssemblyType::BLOCK_SCHUR) {
1214
1219
1220#ifdef ADD_CONTACT
1221
1222 auto get_nested_mat_data = [&](auto schur_dm, auto block_dm) {
1225
1226 {
1227
1228 {simple->getDomainFEName(),
1229
1230 {{"U", "U"},
1231 {"SIGMA", "SIGMA"},
1232 {"U", "SIGMA"},
1233 {"SIGMA", "U"},
1234 {"EP", "EP"},
1235 {"TAU", "TAU"},
1236 {"U", "EP"},
1237 {"EP", "U"},
1238 {"EP", "TAU"},
1239 {"TAU", "EP"},
1240 {"TAU", "U"}
1241
1242 }},
1243
1244 {simple->getBoundaryFEName(),
1245
1246 {{"SIGMA", "SIGMA"}, {"U", "SIGMA"}, {"SIGMA", "U"}
1247
1248 }}
1249
1250 }
1251
1252 );
1253
1255
1256 {dm_schur, dm_block}, block_mat_data,
1257
1258 {"SIGMA", "EP", "TAU"}, {nullptr, nullptr, nullptr}, true
1259
1260 );
1261 };
1262
1263#else
1264
1265 auto get_nested_mat_data = [&](auto schur_dm, auto block_dm) {
1266 auto block_mat_data =
1268
1269 {{simple->getDomainFEName(),
1270
1271 {{"U", "U"},
1272 {"EP", "EP"},
1273 {"TAU", "TAU"},
1274 {"U", "EP"},
1275 {"EP", "U"},
1276 {"EP", "TAU"},
1277 {"TAU", "U"},
1278 {"TAU", "EP"}
1279
1280 }}}
1281
1282 );
1283
1285
1286 {dm_schur, dm_block}, block_mat_data,
1287
1288 {"EP", "TAU"}, {nullptr, nullptr}, false
1289
1290 );
1291 };
1292
1293#endif
1294
1295 auto nested_mat_data = get_nested_mat_data(dm_schur, dm_block);
1297
1298 auto block_is =
getDMSubData(dm_block)->getSmartRowIs();
1299 auto ao_schur =
getDMSubData(dm_schur)->getSmartRowMap();
1300
1301
1302
1303 schur_ptr =
1305 CHKERR schur_ptr->setUp(solver);
1306 }
1307
1309 };
1310
1314 CHKERR VecSetDM(
D, PETSC_NULLPTR);
1315 CHKERR VecSetDM(DD, PETSC_NULLPTR);
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>();
1407 ts_ctx_ptr->getPreProcessIFunction().push_front(pre_proc_ptr);
1408 ts_ctx_ptr->getPreProcessIJacobian().push_front(pre_proc_ptr);
1409 ts_ctx_ptr->getPostProcessIFunction().push_back(post_proc_rhs_ptr);
1410 ts_ctx_ptr->getPostProcessIJacobian().push_back(post_proc_lhs_ptr);
1411
1412
1413 auto disp_time_scale = boost::make_shared<TimeScale>();
1414
1415 auto get_bc_hook_rhs = [&]() {
1417 mField, pre_proc_ptr, {disp_time_scale},
false);
1418 };
1419 pre_proc_ptr->preProcessHook = get_bc_hook_rhs();
1420
1421 auto waak_post_proc_rhs_ptr = boost::weak_ptr<FEMethod>(
1422 post_proc_rhs_ptr);
1423 auto get_post_proc_hook_rhs = [this, waak_post_proc_rhs_ptr]() {
1426 mField, waak_post_proc_rhs_ptr.lock(),
nullptr, Sev::verbose)();
1428 mField, waak_post_proc_rhs_ptr.lock(), 1.)();
1430 };
1431 auto get_post_proc_hook_lhs = [&]() {
1433 mField, post_proc_lhs_ptr, 1.);
1434 };
1435
1436 post_proc_rhs_ptr->postProcessHook = get_post_proc_hook_rhs;
1437 post_proc_lhs_ptr->postProcessHook = get_post_proc_hook_lhs();
1438
1440 };
1441
1444 CHKERR TSSetIJacobian(solver,
B,
B, PETSC_NULLPTR, PETSC_NULLPTR);
1445 } else {
1446 CHKERR TSSetI2Jacobian(solver,
B,
B, PETSC_NULLPTR, PETSC_NULLPTR);
1447 }
1449 CHKERR TSSetSolution(solver,
D);
1450 } else {
1451 CHKERR TS2SetSolution(solver,
D, DD);
1452 }
1453 CHKERR set_section_monitor(solver);
1454 CHKERR set_time_monitor(dm, solver);
1455 CHKERR TSSetFromOptions(solver);
1456
1457 CHKERR add_active_dofs_elem(dm);
1458 boost::shared_ptr<SetUpSchur> schur_ptr;
1459 CHKERR set_schur_pc(solver, schur_ptr);
1460 CHKERR set_essential_bc(dm, solver);
1461
1465 BOOST_LOG_SCOPED_THREAD_ATTR("Timeline", attrs::timer());
1466 MOFEM_LOG(
"TIMER", Sev::verbose) <<
"TSSetUp";
1468 MOFEM_LOG(
"TIMER", Sev::verbose) <<
"TSSetUp <= done";
1469 MOFEM_LOG(
"TIMER", Sev::verbose) <<
"TSSolve";
1470 CHKERR TSSolve(solver, NULL);
1471 MOFEM_LOG(
"TIMER", Sev::verbose) <<
"TSSolve <= done";
1472
1476 "ts_manager_graph.dot");
1477 }
1478
1480}
#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)
