795 {
799
800 auto dm =
simple->getDM();
801 auto time_scale = boost::make_shared<TimeScale>();
803
804
805 auto create_post_proc_fe = [dm,
this,
simple]() {
806
807
808
809 auto post_proc_ele_domain = [this](auto &pip_domain, auto &fe_name) {
810
812 pip_domain, {
H1,
HDIV},
"GEOMETRY");
813
814
815
816
817 auto grad_ptr = boost::make_shared<MatrixDouble>();
818 pip_domain.push_back(
820 grad_ptr));
821
822
823
824
827
828 pip_domain.push_back(op_this);
829
830
831
832 auto fe_physics = op_this->getThisFEPtr();
833
834 auto evaluate_stress_on_physical_element = [&]() {
835
836 fe_physics->getRuleHook =
cpPtr->integrationRule;
838 fe_physics->getOpPtrVector(), {H1});
839 auto common_data_ptr =
840 boost::make_shared<ADOLCPlasticity::CommonData>();
841
842
846 mField,
"U", fe_physics->getOpPtrVector(),
"ADOLCMAT", common_data_ptr,
cpPtr);
847
848 } else {
849 fe_physics->getOpPtrVector().push_back(
851 "U", common_data_ptr->getGradAtGaussPtsPtr()));
852 CHKERR cpPtr->addMatBlockOps(
mField, fe_physics->getOpPtrVector(),
"ADOLCMAT", Sev::inform);
855 }
856 return common_data_ptr;
857 };
858
859 auto dg_projection_froward_and_back = [&](auto &common_data_ptr) {
860
862 auto entity_data_l2 =
863 boost::make_shared<EntitiesFieldData>(MBENTITYSET);
864 auto mass_ptr =
865 boost::make_shared<MatrixDouble>();
866
869
870 auto coeffs_ptr_stress = boost::make_shared<MatrixDouble>();
871
873 common_data_ptr->getStressMatrixPtr(), coeffs_ptr_stress, mass_ptr,
875
876 auto coeffs_ptr_plastic_strain = boost::make_shared<MatrixDouble>();
878 common_data_ptr->getPlasticStrainMatrixPtr(),
879 coeffs_ptr_plastic_strain, mass_ptr, entity_data_l2,
approxBase,
881
882
883
885 common_data_ptr->getStressMatrixPtr(), coeffs_ptr_stress,
888 common_data_ptr->getPlasticStrainMatrixPtr(),
889 coeffs_ptr_plastic_strain, entity_data_l2,
approxBase,
L2));
890 };
891
892 auto common_data_ptr = evaluate_stress_on_physical_element();
893 dg_projection_froward_and_back(common_data_ptr);
894
895 return boost::make_tuple(grad_ptr, common_data_ptr->getStressMatrixPtr(),
896 common_data_ptr->getPlasticStrainMatrixPtr());
897 };
898
899
900
901
902 auto add_post_proc_map = [&](auto post_proc_fe, auto u_ptr, auto grad_ptr,
903 auto stress_ptr, auto plastic_strain_ptr,
904 auto contact_stress_ptr, auto X_ptr) {
906 post_proc_fe->getOpPtrVector().push_back(
907
908 new OpPPMapSPACE_DIM(
909
910 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
911
912 {},
913
914 {{"U", u_ptr}, {"GEOMETRY", X_ptr}},
915
916 {{"GRAD", grad_ptr}, {"SIGMA", contact_stress_ptr}},
917
918 {}
919
920 )
921
922 );
923
926 post_proc_fe->getOpPtrVector().push_back(
927
928 new OpPPMap3D(
929
930 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
931
932 {},
933
934 {},
935
936 {{"PK1", stress_ptr}},
937
938 {{"PLASTIC_STRAIN", plastic_strain_ptr}}
939
940 )
941
942 );
943 } else {
944 post_proc_fe->getOpPtrVector().push_back(
945
946 new OpPPMap3D(
947
948 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
949
950 {},
951
952 {},
953
954 {},
955
956 {{"STRESS", stress_ptr}, {"PLASTIC_STRAIN", plastic_strain_ptr}}
957
958 )
959
960 );
961 }
962
963 return post_proc_fe;
964 };
965
966 auto vol_post_proc = [
this,
simple, post_proc_ele_domain,
967 add_post_proc_map]() {
968 PetscBool post_proc_vol = PETSC_FALSE;
970 post_proc_vol = PETSC_TRUE;
971
973 &post_proc_vol, PETSC_NULLPTR);
974 if (post_proc_vol == PETSC_FALSE)
975 return boost::shared_ptr<PostProcEleDomain>();
976 auto post_proc_fe = boost::make_shared<PostProcEleDomain>(
mField);
977 auto u_ptr = boost::make_shared<MatrixDouble>();
978 post_proc_fe->getOpPtrVector().push_back(
980 auto X_ptr = boost::make_shared<MatrixDouble>();
981 post_proc_fe->getOpPtrVector().push_back(
983 auto contact_stress_ptr = boost::make_shared<MatrixDouble>();
984#ifdef ADD_CONTACT
985 post_proc_fe->getOpPtrVector().push_back(
987 "SIGMA", contact_stress_ptr));
988#else
989 contact_stress_ptr = nullptr;
990#endif
991 auto [grad_ptr, stress_ptr, plastic_strain_ptr] = post_proc_ele_domain(
992 post_proc_fe->getOpPtrVector(),
simple->getDomainFEName());
993
994 return add_post_proc_map(post_proc_fe, u_ptr, grad_ptr, stress_ptr,
995 plastic_strain_ptr, contact_stress_ptr, X_ptr);
996 };
997
998 auto skin_post_proc = [
this,
simple, post_proc_ele_domain,
999 add_post_proc_map]() {
1000
1001
1002 PetscBool post_proc_skin = PETSC_TRUE;
1004 post_proc_skin = PETSC_FALSE;
1005
1007 &post_proc_skin, PETSC_NULLPTR);
1008
1009 if (post_proc_skin == PETSC_FALSE)
1010 return boost::shared_ptr<PostProcEleBdy>();
1011
1012 auto post_proc_fe = boost::make_shared<PostProcEleBdy>(
mField);
1013 auto u_ptr = boost::make_shared<MatrixDouble>();
1014 post_proc_fe->getOpPtrVector().push_back(
1016 auto X_ptr = boost::make_shared<MatrixDouble>();
1017 post_proc_fe->getOpPtrVector().push_back(
1019 auto contact_stress_ptr = boost::make_shared<MatrixDouble>();
1020 auto op_loop_side =
1022 auto [grad_ptr, stress_ptr, plastic_strain_ptr] = post_proc_ele_domain(
1023 op_loop_side->getOpPtrVector(),
simple->getDomainFEName());
1024#ifdef ADD_CONTACT
1025 op_loop_side->getOpPtrVector().push_back(
1027 "SIGMA", contact_stress_ptr));
1028#else
1029 contact_stress_ptr = nullptr;
1030#endif
1031 post_proc_fe->getOpPtrVector().push_back(op_loop_side);
1032
1033 return add_post_proc_map(post_proc_fe, u_ptr, grad_ptr, stress_ptr,
1034 plastic_strain_ptr, contact_stress_ptr, X_ptr);
1035 };
1036
1037 return std::make_pair(vol_post_proc(), skin_post_proc());
1038 };
1039
1040 auto create_reaction_fe = [&]() {
1041 auto fe_ptr = boost::make_shared<DomainEle>(
mField);
1042 fe_ptr->getRuleHook =
cpPtr->integrationRule;
1043
1044 auto &pip = fe_ptr->getOpPtrVector();
1046 auto common_data_ptr = boost::make_shared<ADOLCPlasticity::CommonData>();
1047
1051 mField,
"U", pip,
"ADOLCMAT", common_data_ptr,
cpPtr);
1056 pip.push_back(
1058 } else {
1060 "U", common_data_ptr->getGradAtGaussPtsPtr()));
1062 mField,
"U", pip,
"ADOLCMAT", common_data_ptr,
cpPtr);
1063 }
1064
1065 return fe_ptr;
1066 };
1067
1068 auto add_extra_finite_elements_to_ksp_solver_pipelines = [&]() {
1070
1071 auto pre_proc_ptr = boost::make_shared<FEMethod>();
1072 auto post_proc_rhs_ptr = boost::make_shared<FEMethod>();
1073 auto post_proc_lhs_ptr = boost::make_shared<FEMethod>();
1074
1075 auto get_bc_hook_rhs = [this, pre_proc_ptr, time_scale]() {
1078 {time_scale}, false)();
1080 };
1081
1082 pre_proc_ptr->preProcessHook = get_bc_hook_rhs;
1083
1084 auto get_post_proc_hook_rhs = [this, post_proc_rhs_ptr]() {
1087 mField, post_proc_rhs_ptr,
nullptr, Sev::verbose)();
1089 mField, post_proc_rhs_ptr, 1.)();
1091 };
1092 auto get_post_proc_hook_lhs = [this, post_proc_lhs_ptr]() {
1095 mField, post_proc_lhs_ptr, 1.)();
1097 };
1098 post_proc_rhs_ptr->postProcessHook = get_post_proc_hook_rhs;
1099 post_proc_lhs_ptr->postProcessHook = get_post_proc_hook_lhs;
1100
1101
1103 ts_ctx_ptr->getPreProcessIFunction().push_front(pre_proc_ptr);
1104 ts_ctx_ptr->getPreProcessIJacobian().push_front(pre_proc_ptr);
1105 ts_ctx_ptr->getPostProcessIFunction().push_back(post_proc_rhs_ptr);
1106 ts_ctx_ptr->getPostProcessIJacobian().push_back(post_proc_lhs_ptr);
1108 };
1109
1110
1111
1112 CHKERR add_extra_finite_elements_to_ksp_solver_pipelines();
1113
1114 auto create_monitor_fe = [dm, time_scale](auto &&post_proc_fe,
1115 auto &&reaction_fe) {
1116 return boost::make_shared<Monitor>(
1117 dm, post_proc_fe, reaction_fe,
1118 std::vector<boost::shared_ptr<ScalingMethod>>{time_scale});
1119 };
1120
1121
1122 auto set_up_post_step = [&](auto ts) {
1124
1125
1126
1127 auto create_update_ptr = [&]() {
1128 auto fe_ptr = boost::make_shared<DomainEle>(
mField);
1129 fe_ptr->getRuleHook =
cpPtr->integrationRule;
1131 {H1});
1132 auto common_data_ptr = boost::make_shared<ADOLCPlasticity::CommonData>();
1133
1137 mField,
"U", fe_ptr->getOpPtrVector(),
"ADOLCMAT", common_data_ptr,
1139 } else {
1140 fe_ptr->getOpPtrVector().push_back(
1142 "U", common_data_ptr->getGradAtGaussPtsPtr()));
1144 "ADOLCMAT", Sev::noisy);
1145 fe_ptr->getOpPtrVector().push_back(
1148 }
1149
1152 "ADOLCMAT", common_data_ptr,
cpPtr),
1153 "push update ops");
1154 return fe_ptr;
1155 };
1156
1157
1160
1161
1162
1163 auto ts_step_post_proc = [](TS ts) {
1168 };
1169
1170
1171
1172 CHKERR TSSetPostStep(ts, ts_step_post_proc);
1174 };
1175
1176
1177
1178 auto set_up_monitor = [&](auto ts) {
1180 boost::shared_ptr<FEMethod> null_fe;
1181 auto monitor_ptr =
1182 create_monitor_fe(create_post_proc_fe(), create_reaction_fe());
1184 null_fe, monitor_ptr);
1186 };
1187
1188 auto set_section_monitor = [&](auto solver) {
1190 SNES snes;
1191 CHKERR TSGetSNES(solver, &snes);
1192 CHKERR SNESMonitorSet(snes,
1195 (void *)(snes_ctx_ptr.get()), nullptr);
1197 };
1198
1199 auto set_up_adapt = [&](auto ts) {
1202 TSAdapt adapt;
1203 CHKERR TSGetAdapt(ts, &adapt);
1205 };
1206
1207
1208 auto ts = pipeline_mng->createTSIM();
1209
1210
1211 double ftime = 1;
1212 CHKERR TSSetMaxTime(ts, ftime);
1213 CHKERR TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);
1216 CHKERR TSSetIJacobian(ts,
B,
B, PETSC_NULLPTR, PETSC_NULLPTR);
1217#ifdef ADD_CONTACT
1219#endif
1221 CHKERR set_section_monitor(ts);
1222
1223
1224 CHKERR set_up_monitor(ts);
1225 CHKERR set_up_post_step(ts);
1227 CHKERR TSSetFromOptions(ts);
1228
1229 CHKERR TSSolve(ts, NULL);
1230
1231
1232 CHKERR TSGetTime(ts, &ftime);
1233
1234 PetscInt steps, snesfails, rejects, nonlinits, linits;
1235 CHKERR TSGetStepNumber(ts, &steps);
1236 CHKERR TSGetSNESFailures(ts, &snesfails);
1237 CHKERR TSGetStepRejections(ts, &rejects);
1238 CHKERR TSGetSNESIterations(ts, &nonlinits);
1239 CHKERR TSGetKSPIterations(ts, &linits);
1241 "steps %d (%d rejected, %d SNES fails), ftime %g, nonlinits "
1242 "%d, linits %d",
1243 steps, rejects, snesfails, ftime, nonlinits, linits);
1244
1246}
#define MOFEM_LOG_C(channel, severity, format,...)
static boost::shared_ptr< TSUpdate > ts_update_ptr
PipelineManager::ElementsAndOpsByDim< SPACE_DIM >::DomainEle DomainEle
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
@ L2
field with C-1 continuity
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
MoFEMErrorCode opFactoryDomainUpdate(MoFEM::Interface &m_field, Pip &pip, std::string block_name, boost::shared_ptr< CommonData > common_data_ptr, boost::shared_ptr< ClosestPointProjection > cp_ptr)
Push operators to update history variables.
auto createDMMatrix(DM dm)
Get smart matrix from DM.
ForcesAndSourcesCore::UserDataOperator * getRawPtrOpCalculateStress(MoFEM::Interface &m_field, boost::shared_ptr< CommonData > common_data_ptr, boost::shared_ptr< ClosestPointProjection > cp_ptr, bool calc_lhs)
Get opreator to calulate stress.
boost::shared_ptr< TSUpdate > createTSUpdate(std::string fe_name, boost::shared_ptr< FEMethod > fe_ptr)
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
PetscErrorCode DMMoFEMTSSetMonitor(DM dm, TS ts, const std::string fe_name, boost::shared_ptr< MoFEM::FEMethod > method, boost::shared_ptr< MoFEM::BasicMethod > pre_only, boost::shared_ptr< MoFEM::BasicMethod > post_only)
Set Monitor To TS solver.
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 getDMSnesCtx(DM dm)
Get SNES context data structure used by DM.
PetscErrorCode TSAdaptCreateMoFEM(TSAdapt adapt)
Craete MOFEM adapt.
Class (Function) to enforce essential constrains on the left hand side diagonal.
Class (Function) to enforce essential constrains on the right hand side diagonal.
Class (Function) to calculate residual side diagonal.
Class (Function) to enforce essential constrains.
Calculate tenor field using vectorial base, i.e. Hdiv/Hcurl.
Element used to execute operators on side of the element.
Execute "this" element in the operator.
Post post-proc data at points from hash maps.
