827 auto set_section_monitor = [&](
auto solver) {
830 CHKERR TSGetSNES(solver, &snes);
831 CHKERR SNESMonitorSet(snes,
834 (
void *)(snes_ctx_ptr.get()),
nullptr);
838 auto create_post_process_elements = [&]() {
839 auto push_vol_ops = [
this](
auto &pip) {
841 pip, {
H1,
HDIV},
"GEOMETRY");
843 auto [common_plastic_ptr, common_hencky_ptr] =
844 PlasticOps::createCommonPlasticOps<SPACE_DIM, IT, DomainEleOp>(
845 mField,
"MAT_PLASTIC", pip,
"U",
"EP",
"TAU", 1., Sev::inform);
847 if (common_hencky_ptr) {
848 if (common_plastic_ptr->mGradPtr != common_hencky_ptr->matGradPtr)
852 return std::make_pair(common_plastic_ptr, common_hencky_ptr);
855 auto push_vol_post_proc_ops = [
this](
auto &pp_fe,
auto &&p) {
858 auto &pip = pp_fe->getOpPtrVector();
860 auto [common_plastic_ptr, common_hencky_ptr] = p;
864 auto x_ptr = boost::make_shared<MatrixDouble>();
867 auto u_ptr = boost::make_shared<MatrixDouble>();
876 pp_fe->getPostProcMesh(), pp_fe->getMapGaussPts(),
879 common_plastic_ptr->getPlasticSurfacePtr()},
880 {
"PLASTIC_MULTIPLIER",
881 common_plastic_ptr->getPlasticTauPtr()}},
883 {{
"U", u_ptr}, {
"GEOMETRY", x_ptr}},
885 {{
"GRAD", common_hencky_ptr->matGradPtr},
886 {
"FIRST_PIOLA", common_hencky_ptr->getMatFirstPiolaStress()}},
888 {{
"HENCKY_STRAIN", common_hencky_ptr->getMatLogC()},
889 {
"PLASTIC_STRAIN", common_plastic_ptr->getPlasticStrainPtr()},
890 {
"PLASTIC_FLOW", common_plastic_ptr->getPlasticFlowPtr()}}
902 pp_fe->getPostProcMesh(), pp_fe->getMapGaussPts(),
905 common_plastic_ptr->getPlasticSurfacePtr()},
906 {
"PLASTIC_MULTIPLIER",
907 common_plastic_ptr->getPlasticTauPtr()}},
909 {{
"U", u_ptr}, {
"GEOMETRY", x_ptr}},
913 {{
"STRAIN", common_plastic_ptr->mStrainPtr},
914 {
"STRESS", common_plastic_ptr->mStressPtr},
915 {
"PLASTIC_STRAIN", common_plastic_ptr->getPlasticStrainPtr()},
916 {
"PLASTIC_FLOW", common_plastic_ptr->getPlasticFlowPtr()}}
926 PetscBool post_proc_vol;
927 PetscBool post_proc_skin;
930 post_proc_vol = PETSC_TRUE;
931 post_proc_skin = PETSC_FALSE;
933 post_proc_vol = PETSC_FALSE;
934 post_proc_skin = PETSC_TRUE;
939 &post_proc_skin, PETSC_NULL);
941 auto vol_post_proc = [
this, push_vol_post_proc_ops, push_vol_ops,
943 if (post_proc_vol == PETSC_FALSE)
944 return boost::shared_ptr<PostProcEle>();
945 auto pp_fe = boost::make_shared<PostProcEle>(
mField);
947 push_vol_post_proc_ops(pp_fe, push_vol_ops(pp_fe->getOpPtrVector())),
948 "push_vol_post_proc_ops");
952 auto skin_post_proc = [
this, push_vol_post_proc_ops, push_vol_ops,
954 if (post_proc_skin == PETSC_FALSE)
955 return boost::shared_ptr<SkinPostProcEle>();
958 auto pp_fe = boost::make_shared<SkinPostProcEle>(
mField);
961 pp_fe->getOpPtrVector().push_back(op_side);
963 pp_fe, push_vol_ops(op_side->getOpPtrVector())),
964 "push_vol_post_proc_ops");
968 return std::make_pair(vol_post_proc(), skin_post_proc());
971 auto scatter_create = [&](
auto D,
auto coeff) {
973 CHKERR is_manager->isCreateProblemFieldAndRank(
simple->getProblemName(),
974 ROW,
"U", coeff, coeff, is);
976 CHKERR ISGetLocalSize(is, &loc_size);
980 CHKERR VecScatterCreate(
D, is,
v, PETSC_NULL, &scatter);
985 boost::shared_ptr<SetPtsData> field_eval_data;
986 boost::shared_ptr<MatrixDouble> u_field_ptr;
988 std::array<double, SPACE_DIM> field_eval_coords;
991 field_eval_coords.data(), &coords_dim,
994 boost::shared_ptr<std::map<std::string, boost::shared_ptr<VectorDouble>>>
995 scalar_field_ptrs = boost::make_shared<
996 std::map<std::string, boost::shared_ptr<VectorDouble>>>();
997 boost::shared_ptr<std::map<std::string, boost::shared_ptr<MatrixDouble>>>
998 vector_field_ptrs = boost::make_shared<
999 std::map<std::string, boost::shared_ptr<MatrixDouble>>>();
1000 boost::shared_ptr<std::map<std::string, boost::shared_ptr<MatrixDouble>>>
1001 sym_tensor_field_ptrs = boost::make_shared<
1002 std::map<std::string, boost::shared_ptr<MatrixDouble>>>();
1003 boost::shared_ptr<std::map<std::string, boost::shared_ptr<MatrixDouble>>>
1004 tensor_field_ptrs = boost::make_shared<
1005 std::map<std::string, boost::shared_ptr<MatrixDouble>>>();
1008 auto u_field_ptr = boost::make_shared<MatrixDouble>();
1013 field_eval_data,
simple->getDomainFEName());
1015 field_eval_data->setEvalPoints(field_eval_coords.data(), 1);
1016 auto no_rule = [](
int,
int,
int) {
return -1; };
1017 auto field_eval_fe_ptr = field_eval_data->feMethodPtr.lock();
1018 field_eval_fe_ptr->getRuleHook = no_rule;
1021 field_eval_fe_ptr->getOpPtrVector(), {H1, HDIV},
"GEOMETRY");
1023 auto [common_plastic_ptr, common_hencky_ptr] =
1024 PlasticOps::createCommonPlasticOps<SPACE_DIM, IT, DomainEleOp>(
1025 mField,
"MAT_PLASTIC", field_eval_fe_ptr->getOpPtrVector(),
"U",
1026 "EP",
"TAU", 1., Sev::inform);
1028 field_eval_fe_ptr->getOpPtrVector().push_back(
1031 if ((common_plastic_ptr) && (common_hencky_ptr) && (scalar_field_ptrs)) {
1033 scalar_field_ptrs->insert(
1034 {
"PLASTIC_SURFACE", common_plastic_ptr->getPlasticSurfacePtr()});
1035 scalar_field_ptrs->insert(
1036 {
"PLASTIC_MULTIPLIER", common_plastic_ptr->getPlasticTauPtr()});
1037 vector_field_ptrs->insert({
"U", u_field_ptr});
1038 sym_tensor_field_ptrs->insert(
1039 {
"PLASTIC_STRAIN", common_plastic_ptr->getPlasticStrainPtr()});
1040 sym_tensor_field_ptrs->insert(
1041 {
"PLASTIC_FLOW", common_plastic_ptr->getPlasticFlowPtr()});
1042 sym_tensor_field_ptrs->insert(
1043 {
"HENCKY_STRAIN", common_hencky_ptr->getMatLogC()});
1044 tensor_field_ptrs->insert({
"GRAD", common_hencky_ptr->matGradPtr});
1045 tensor_field_ptrs->insert(
1046 {
"FIRST_PIOLA", common_hencky_ptr->getMatFirstPiolaStress()});
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 {
"STRAIN", common_plastic_ptr->mStrainPtr});
1055 sym_tensor_field_ptrs->insert(
1056 {
"STRESS", common_plastic_ptr->mStressPtr});
1057 sym_tensor_field_ptrs->insert(
1058 {
"PLASTIC_STRAIN", common_plastic_ptr->getPlasticStrainPtr()});
1059 sym_tensor_field_ptrs->insert(
1060 {
"PLASTIC_FLOW", common_plastic_ptr->getPlasticFlowPtr()});
1065 auto test_monitor_ptr = boost::make_shared<FEMethod>();
1067 auto set_time_monitor = [&](
auto dm,
auto solver) {
1071 uZScatter, field_eval_coords, field_eval_data, scalar_field_ptrs,
1072 vector_field_ptrs, sym_tensor_field_ptrs, tensor_field_ptrs));
1073 boost::shared_ptr<ForcesAndSourcesCore>
null;
1075 test_monitor_ptr->postProcessHook = [&]() {
1078 if (
atom_test && fabs(test_monitor_ptr->ts_t - 0.5) < 1e-12 &&
1079 test_monitor_ptr->ts_step == 25) {
1081 if (scalar_field_ptrs->at(
"PLASTIC_MULTIPLIER")->size()) {
1084 MOFEM_LOG(
"PlasticSync", Sev::inform) <<
"Eval point tau: " << t_tau;
1086 if (
atom_test == 1 && fabs(t_tau - 0.688861) > 1e-5) {
1088 "atom test %d failed: wrong plastic multiplier value",
1093 if (vector_field_ptrs->at(
"U")->size1()) {
1096 getFTensor1FromMat<SPACE_DIM>(*vector_field_ptrs->at(
"U"));
1097 MOFEM_LOG(
"PlasticSync", Sev::inform) <<
"Eval point U: " << t_disp;
1099 if (
atom_test == 1 && fabs(t_disp(0) - 0.25 / 2.) > 1e-5 ||
1100 fabs(t_disp(1) + 0.0526736) > 1e-5 || fabs(t_disp(2)) > 1e-5) {
1102 "atom test %d failed: wrong displacement value",
1107 if (sym_tensor_field_ptrs->at(
"PLASTIC_STRAIN")->size1()) {
1108 auto t_plastic_strain = getFTensor2SymmetricFromMat<SPACE_DIM>(
1109 *sym_tensor_field_ptrs->at(
"PLASTIC_STRAIN"));
1111 <<
"Eval point EP: " << t_plastic_strain;
1114 fabs(t_plastic_strain(0, 0) - 0.221943) > 1e-5 ||
1115 fabs(t_plastic_strain(0, 1)) > 1e-5 ||
1116 fabs(t_plastic_strain(1, 1) + 0.110971) > 1e-5) {
1118 "atom test %d failed: wrong plastic strain value",
1123 if (tensor_field_ptrs->at(
"FIRST_PIOLA")->size1()) {
1124 auto t_piola_stress = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(
1125 *tensor_field_ptrs->at(
"FIRST_PIOLA"));
1127 <<
"Eval point Piola stress: " << t_piola_stress;
1129 if (
atom_test == 1 && fabs((t_piola_stress(0, 0) - 198.775) /
1130 t_piola_stress(0, 0)) > 1e-5 ||
1131 fabs(t_piola_stress(0, 1)) + fabs(t_piola_stress(1, 0)) +
1132 fabs(t_piola_stress(1, 1)) >
1135 "atom test %d failed: wrong Piola stress value",
1146 monitor_ptr,
null, test_monitor_ptr);
1151 auto set_schur_pc = [&](
auto solver,
1152 boost::shared_ptr<SetUpSchur> &schur_ptr) {
1155 auto name_prb =
simple->getProblemName();
1166 for (
auto f : {
"U"}) {
1184 for (
auto f : {
"SIGMA",
"EP",
"TAU"}) {
1189 for (
auto f : {
"EP",
"TAU"}) {
1208 auto get_nested_mat_data = [&](
auto schur_dm,
auto block_dm) {
1214 {simple->getDomainFEName(),
1230 {simple->getBoundaryFEName(),
1232 {{
"SIGMA",
"SIGMA"}, {
"U",
"SIGMA"}, {
"SIGMA",
"U"}
1242 {dm_schur, dm_block}, block_mat_data,
1244 {
"SIGMA",
"EP",
"TAU"}, {
nullptr,
nullptr,
nullptr}, true
1251 auto get_nested_mat_data = [&](
auto schur_dm,
auto block_dm) {
1252 auto block_mat_data =
1255 {{simple->getDomainFEName(),
1272 {dm_schur, dm_block}, block_mat_data,
1274 {
"EP",
"TAU"}, {
nullptr,
nullptr}, false
1281 auto nested_mat_data = get_nested_mat_data(dm_schur, dm_block);
1284 auto block_is =
getDMSubData(dm_block)->getSmartRowIs();
1285 auto ao_schur =
getDMSubData(dm_schur)->getSmartRowMap();
1291 CHKERR schur_ptr->setUp(solver);
1297 auto dm =
simple->getDM();
1305 auto create_solver = [pip_mng]() {
1307 return pip_mng->createTSIM();
1309 return pip_mng->createTSIM2();
1312 auto solver = create_solver();
1314 auto active_pre_lhs = []() {
1321 auto active_post_lhs = [&]() {
1323 auto get_iter = [&]() {
1328 "Can not get iter");
1332 auto iter = get_iter();
1335 std::array<int, 5> activity_data;
1336 std::fill(activity_data.begin(), activity_data.end(), 0);
1338 activity_data.data(), activity_data.size(), MPI_INT,
1341 int &active_points = activity_data[0];
1342 int &avtive_full_elems = activity_data[1];
1343 int &avtive_elems = activity_data[2];
1344 int &nb_points = activity_data[3];
1345 int &nb_elements = activity_data[4];
1349 double proc_nb_points =
1350 100 *
static_cast<double>(active_points) / nb_points;
1351 double proc_nb_active =
1352 100 *
static_cast<double>(avtive_elems) / nb_elements;
1353 double proc_nb_full_active = 100;
1355 proc_nb_full_active =
1356 100 *
static_cast<double>(avtive_full_elems) / avtive_elems;
1359 "Iter %d nb pts %d nb active pts %d (%3.3f\%) nb active "
1361 "(%3.3f\%) nb full active elems %d (%3.3f\%)",
1362 iter, nb_points, active_points, proc_nb_points,
1363 avtive_elems, proc_nb_active, avtive_full_elems,
1364 proc_nb_full_active, iter);
1371 auto add_active_dofs_elem = [&](
auto dm) {
1373 auto fe_pre_proc = boost::make_shared<FEMethod>();
1374 fe_pre_proc->preProcessHook = active_pre_lhs;
1375 auto fe_post_proc = boost::make_shared<FEMethod>();
1376 fe_post_proc->postProcessHook = active_post_lhs;
1378 ts_ctx_ptr->getPreProcessIJacobian().push_front(fe_pre_proc);
1379 ts_ctx_ptr->getPostProcessIJacobian().push_back(fe_post_proc);
1383 auto set_essential_bc = [&](
auto dm,
auto solver) {
1387 auto pre_proc_ptr = boost::make_shared<FEMethod>();
1388 auto post_proc_rhs_ptr = boost::make_shared<FEMethod>();
1389 auto post_proc_lhs_ptr = boost::make_shared<FEMethod>();
1392 auto disp_time_scale = boost::make_shared<TimeScale>();
1394 auto get_bc_hook_rhs = [
this, pre_proc_ptr, disp_time_scale]() {
1396 {disp_time_scale},
false);
1399 pre_proc_ptr->preProcessHook = get_bc_hook_rhs();
1401 auto get_post_proc_hook_rhs = [
this, post_proc_rhs_ptr]() {
1404 mField, post_proc_rhs_ptr,
nullptr, Sev::verbose)();
1406 mField, post_proc_rhs_ptr, 1.)();
1409 auto get_post_proc_hook_lhs = [
this, post_proc_lhs_ptr]() {
1411 mField, post_proc_lhs_ptr, 1.);
1413 post_proc_rhs_ptr->postProcessHook = get_post_proc_hook_rhs;
1416 ts_ctx_ptr->getPreProcessIFunction().push_front(pre_proc_ptr);
1417 ts_ctx_ptr->getPreProcessIJacobian().push_front(pre_proc_ptr);
1418 ts_ctx_ptr->getPostProcessIFunction().push_back(post_proc_rhs_ptr);
1419 post_proc_lhs_ptr->postProcessHook = get_post_proc_hook_lhs();
1420 ts_ctx_ptr->getPostProcessIJacobian().push_back(post_proc_lhs_ptr);
1426 CHKERR TSSetSolution(solver,
D);
1428 CHKERR TS2SetSolution(solver,
D, DD);
1431 CHKERR set_section_monitor(solver);
1432 CHKERR set_time_monitor(dm, solver);
1433 CHKERR TSSetFromOptions(solver);
1436 CHKERR add_active_dofs_elem(dm);
1437 boost::shared_ptr<SetUpSchur> schur_ptr;
1438 CHKERR set_schur_pc(solver, schur_ptr);
1439 CHKERR set_essential_bc(dm, solver);
1444 BOOST_LOG_SCOPED_THREAD_ATTR(
"Timeline", attrs::timer());
1445 MOFEM_LOG(
"TIMER", Sev::verbose) <<
"TSSetUp";
1447 MOFEM_LOG(
"TIMER", Sev::verbose) <<
"TSSetUp <= done";
1448 MOFEM_LOG(
"TIMER", Sev::verbose) <<
"TSSolve";
1449 CHKERR TSSolve(solver, NULL);
1450 MOFEM_LOG(
"TIMER", Sev::verbose) <<
"TSSolve <= done";
