778 {
780
784
786
787 auto set_section_monitor = [&](auto solver) {
789 SNES snes;
790 CHKERR TSGetSNES(solver, &snes);
791 CHKERR SNESMonitorSet(snes,
794 (void *)(snes_ctx_ptr.get()), nullptr);
796 };
797
798 auto create_post_process_elements = [&]() {
799 auto pp_fe = boost::make_shared<PostProcEle>(
mField);
800 auto &pip = pp_fe->getOpPtrVector();
801
802 auto push_vol_ops = [this](auto &pip) {
804 "GEOMETRY");
805
806 auto [common_plastic_ptr, common_henky_ptr] =
807 PlasticOps::createCommonPlasticOps<SPACE_DIM, IT, DomainEleOp>(
808 mField,
"MAT_PLASTIC", pip,
"U",
"EP",
"TAU", 1., Sev::inform);
809
810 if (common_henky_ptr) {
811 if (common_plastic_ptr->mGradPtr != common_henky_ptr->matGradPtr)
813 }
814
815 return std::make_pair(common_plastic_ptr, common_henky_ptr);
816 };
817
818 auto push_vol_post_proc_ops = [
this](
auto &pp_fe,
auto &&
p) {
820
821 auto &pip = pp_fe->getOpPtrVector();
822
823 auto [common_plastic_ptr, common_henky_ptr] =
p;
824
826
827 auto x_ptr = boost::make_shared<MatrixDouble>();
828 pip.push_back(
830 auto u_ptr = boost::make_shared<MatrixDouble>();
832
834
835 pip.push_back(
836
838
839 pp_fe->getPostProcMesh(), pp_fe->getMapGaussPts(),
840
841 {{"PLASTIC_SURFACE",
842 common_plastic_ptr->getPlasticSurfacePtr()},
843 {"PLASTIC_MULTIPLIER",
844 common_plastic_ptr->getPlasticTauPtr()}},
845
846 {{"U", u_ptr}, {"GEOMETRY", x_ptr}},
847
848 {{"GRAD", common_henky_ptr->matGradPtr},
849 {"FIRST_PIOLA", common_henky_ptr->getMatFirstPiolaStress()}},
850
851 {{"PLASTIC_STRAIN", common_plastic_ptr->getPlasticStrainPtr()},
852 {"PLASTIC_FLOW", common_plastic_ptr->getPlasticFlowPtr()}}
853
854 )
855
856 );
857
858 } else {
859
860 pip.push_back(
861
863
864 pp_fe->getPostProcMesh(), pp_fe->getMapGaussPts(),
865
866 {{"PLASTIC_SURFACE",
867 common_plastic_ptr->getPlasticSurfacePtr()},
868 {"PLASTIC_MULTIPLIER",
869 common_plastic_ptr->getPlasticTauPtr()}},
870
871 {{"U", u_ptr}, {"GEOMETRY", x_ptr}},
872
873 {},
874
875 {{"STRAIN", common_plastic_ptr->mStrainPtr},
876 {"STRESS", common_plastic_ptr->mStressPtr},
877 {"PLASTIC_STRAIN", common_plastic_ptr->getPlasticStrainPtr()},
878 {"PLASTIC_FLOW", common_plastic_ptr->getPlasticFlowPtr()}}
879
880 )
881
882 );
883 }
884
886 };
887
888 auto vol_post_proc = [this, push_vol_post_proc_ops, push_vol_ops]() {
889 PetscBool post_proc_vol = PETSC_FALSE;
891 &post_proc_vol, PETSC_NULL);
892 if (post_proc_vol == PETSC_FALSE)
893 return boost::shared_ptr<PostProcEle>();
894 auto pp_fe = boost::make_shared<PostProcEle>(
mField);
896 push_vol_post_proc_ops(pp_fe, push_vol_ops(pp_fe->getOpPtrVector())),
897 "push_vol_post_proc_ops");
898 return pp_fe;
899 };
900
901 auto skin_post_proc = [this, push_vol_post_proc_ops, push_vol_ops]() {
902 PetscBool post_proc_skin = PETSC_TRUE;
904 &post_proc_skin, PETSC_NULL);
905 if (post_proc_skin == PETSC_FALSE)
906 return boost::shared_ptr<SkinPostProcEle>();
907
909 auto pp_fe = boost::make_shared<SkinPostProcEle>(
mField);
912 pp_fe->getOpPtrVector().push_back(op_side);
914 pp_fe, push_vol_ops(op_side->getOpPtrVector())),
915 "push_vol_post_proc_ops");
916 return pp_fe;
917 };
918
919 return std::make_pair(vol_post_proc(), skin_post_proc());
920 };
921
922 auto scatter_create = [&](
auto D,
auto coeff) {
924 CHKERR is_manager->isCreateProblemFieldAndRank(
simple->getProblemName(),
925 ROW,
"U", coeff, coeff, is);
926 int loc_size;
927 CHKERR ISGetLocalSize(is, &loc_size);
930 VecScatter scatter;
931 CHKERR VecScatterCreate(
D, is,
v, PETSC_NULL, &scatter);
934 };
935
936 auto set_time_monitor = [&](auto dm, auto solver) {
938 boost::shared_ptr<Monitor> monitor_ptr(
941 boost::shared_ptr<ForcesAndSourcesCore> null;
943 monitor_ptr, null, null);
945 };
946
947 auto set_schur_pc = [&](auto solver,
948 boost::shared_ptr<SetUpSchur> &schur_ptr) {
950
952 auto name_prb =
simple->getProblemName();
953
954
963 for (
auto f : {
"U"}) {
966 }
968
970 };
971
972
973 if constexpr (
AT == AssemblyType::SCHUR) {
980
981 IS is_union_raw;
982 CHKERR ISExpand(is_epp, is_tau, &is_union_raw);
984
985#ifdef ADD_CONTACT
986 auto add_sigma_to_is = [&](auto is_union) {
988 auto add_sigma_to_is_impl = [&]() {
993 is_sigma);
994 IS is_union_raw_sigma;
995 CHKERR ISExpand(is_union, is_sigma, &is_union_raw_sigma);
998 };
1000 return is_union_sigma;
1001 };
1002 is_union = add_sigma_to_is(is_union);
1003#endif
1004
1007
1008
1009
1012 schur_ptr =
1014 CHKERR schur_ptr->setUp(solver);
1015 }
1016
1018 };
1019
1020 auto dm =
simple->getDM();
1027
1028 auto create_solver = [pip_mng]() {
1030 return pip_mng->createTSIM();
1031 else
1032 return pip_mng->createTSIM2();
1033 };
1034
1035 auto solver = create_solver();
1036
1037 auto active_pre_lhs = []() {
1042 };
1043
1044 auto active_post_lhs = [&]() {
1046 auto get_iter = [&]() {
1047 SNES snes;
1049 int iter;
1051 "Can not get iter");
1052 return iter;
1053 };
1054
1055 auto iter = get_iter();
1056 if (iter >= 0) {
1057
1058 std::array<int, 5> activity_data;
1059 std::fill(activity_data.begin(), activity_data.end(), 0);
1061 activity_data.data(), activity_data.size(), MPI_INT,
1063
1064 int &active_points = activity_data[0];
1065 int &avtive_full_elems = activity_data[1];
1066 int &avtive_elems = activity_data[2];
1067 int &nb_points = activity_data[3];
1068 int &nb_elements = activity_data[4];
1069
1070 if (nb_points) {
1071
1072 double proc_nb_points =
1073 100 * static_cast<double>(active_points) / nb_points;
1074 double proc_nb_active =
1075 100 * static_cast<double>(avtive_elems) / nb_elements;
1076 double proc_nb_full_active = 100;
1077 if (avtive_elems)
1078 proc_nb_full_active =
1079 100 * static_cast<double>(avtive_full_elems) / avtive_elems;
1080
1082 "Iter %d nb pts %d nb avtive pts %d (%3.3f\%) nb active "
1083 "elements %d "
1084 "(%3.3f\%) nb full active elems %d (%3.3f\%)",
1085 iter, nb_points, active_points, proc_nb_points,
1086 avtive_elems, proc_nb_active, avtive_full_elems,
1087 proc_nb_full_active, iter);
1088 }
1089 }
1090
1092 };
1093
1094 auto add_active_dofs_elem = [&](auto dm) {
1096 auto fe_pre_proc = boost::make_shared<FEMethod>();
1097 fe_pre_proc->preProcessHook = active_pre_lhs;
1098 auto fe_post_proc = boost::make_shared<FEMethod>();
1099 fe_post_proc->postProcessHook = active_post_lhs;
1101 ts_ctx_ptr->getPreProcessIJacobian().push_front(fe_pre_proc);
1102 ts_ctx_ptr->getPostProcessIJacobian().push_back(fe_post_proc);
1104 };
1105
1106 auto set_essential_bc = [&](auto dm, auto solver) {
1108
1109
1110 auto pre_proc_ptr = boost::make_shared<FEMethod>();
1111 auto post_proc_rhs_ptr = boost::make_shared<FEMethod>();
1112 auto post_proc_lhs_ptr = boost::make_shared<FEMethod>();
1113
1114
1115 auto disp_time_scale = boost::make_shared<TimeScale>();
1116
1117 auto get_bc_hook_rhs = [this, pre_proc_ptr, disp_time_scale]() {
1119 {disp_time_scale}, false);
1120 return hook;
1121 };
1122 pre_proc_ptr->preProcessHook = get_bc_hook_rhs();
1123
1124 auto get_post_proc_hook_rhs = [this, post_proc_rhs_ptr]() {
1127 mField, post_proc_rhs_ptr,
nullptr, Sev::verbose)();
1129 mField, post_proc_rhs_ptr, 1.)();
1131 };
1132 auto get_post_proc_hook_lhs = [this, post_proc_lhs_ptr]() {
1134 mField, post_proc_lhs_ptr, 1.);
1135 };
1136 post_proc_rhs_ptr->postProcessHook = get_post_proc_hook_rhs;
1137
1139 ts_ctx_ptr->getPreProcessIFunction().push_front(pre_proc_ptr);
1140 ts_ctx_ptr->getPreProcessIJacobian().push_front(pre_proc_ptr);
1141 ts_ctx_ptr->getPostProcessIFunction().push_back(post_proc_rhs_ptr);
1142
1143 SNES snes;
1144 CHKERR TSGetSNES(solver, &snes);
1145 KSP ksp;
1146 CHKERR SNESGetKSP(snes, &ksp);
1147 PC pc;
1148 CHKERR KSPGetPC(ksp, &pc);
1149 PetscBool is_pcfs = PETSC_FALSE;
1150 PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
1151
1152 if (is_pcfs == PETSC_FALSE) {
1153 post_proc_lhs_ptr->postProcessHook = get_post_proc_hook_lhs();
1154 ts_ctx_ptr->getPostProcessIJacobian().push_back(post_proc_lhs_ptr);
1155 }
1157 };
1158
1160 CHKERR TSSetSolution(solver,
D);
1161 } else {
1162 CHKERR TS2SetSolution(solver,
D, DD);
1163 }
1164
1165 CHKERR set_section_monitor(solver);
1166 CHKERR set_time_monitor(dm, solver);
1167 CHKERR TSSetFromOptions(solver);
1169
1170 CHKERR add_active_dofs_elem(dm);
1171 boost::shared_ptr<SetUpSchur> schur_ptr;
1172 CHKERR set_schur_pc(solver, schur_ptr);
1173 CHKERR set_essential_bc(dm, solver);
1174
1178 BOOST_LOG_SCOPED_THREAD_ATTR("Timeline", attrs::timer());
1179 MOFEM_LOG(
"TIMER", Sev::verbose) <<
"TSSetUp";
1181 MOFEM_LOG(
"TIMER", Sev::verbose) <<
"TSSetUp <= done";
1182 MOFEM_LOG(
"TIMER", Sev::verbose) <<
"TSSolve";
1183 CHKERR TSSolve(solver, NULL);
1184 MOFEM_LOG(
"TIMER", Sev::verbose) <<
"TSSolve <= done";
1185
1187}
#define MAX_DOFS_ON_ENTITY
Maximal number of DOFs on entity.
#define CHK_MOAB_THROW(err, msg)
Check error code of MoAB function and throw MoFEM exception.
PetscErrorCode DMMoFEMCreateSubDM(DM subdm, DM dm, const char problem_name[])
Must be called by user to set Sub DM MoFEM data structures.
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
PetscErrorCode DMMoFEMSetSquareProblem(DM dm, PetscBool square_problem)
set squared problem
PetscErrorCode DMMoFEMAddSubFieldRow(DM dm, const char field_name[])
PetscErrorCode DMMoFEMAddSubFieldCol(DM dm, const char field_name[])
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
const FTensor::Tensor2< T, Dim, Dim > Vec
auto getDMTsCtx(DM dm)
Get TS context data structure used by DM.
MoFEMErrorCode MoFEMSNESMonitorFields(SNES snes, PetscInt its, PetscReal fgnorm, SnesCtx *snes_ctx)
Sens monitor printing residual field by field.
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
auto getDMSubData(DM dm)
Get sub problem data structure.
auto createAOMappingIS(IS isapp, IS ispetsc)
Creates an application mapping using two index sets.
auto getDMSnesCtx(DM dm)
Get SNES context data structure used by DM.
auto createDM(MPI_Comm comm, const std::string dm_type_name)
Creates smart DM object.
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uZScatter
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uXScatter
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uYScatter
Class (Function) to enforce essential constrains on the left hand side diagonal.
Class (Function) to enforce essential constrains on the right hand side diagonal.
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, SmartPetscObj< DM > sub_dm, SmartPetscObj< IS > field_split_it, SmartPetscObj< AO > ao_map)
Create data structure for handling Schur complement.
