860 {
862 boost::shared_ptr<TsCtx>
ts_ctx;
864
865 if (std::abs(
alphaRho) > std::numeric_limits<double>::epsilon()) {
866 CHKERR TSSetI2Function(ts,
f, PETSC_NULL, PETSC_NULL);
867 CHKERR TSSetI2Jacobian(ts,
m,
m, PETSC_NULL, PETSC_NULL);
868 } else {
869 CHKERR TSSetIFunction(ts,
f, PETSC_NULL, PETSC_NULL);
870 CHKERR TSSetIJacobian(ts,
m,
m, PETSC_NULL, PETSC_NULL);
871 }
872
874
878 for (auto &fe : list)
879 if (auto fe_cast = dynamic_cast<EpElementBase *>(fe.second))
880 fe_cast->addStreachSchurMatrix(S, aoS);
881 else
884 };
885
889 for (auto &fe : list)
890 if (auto fe_cast = dynamic_cast<EpElementBase *>(fe.getSharedPtr().get()))
891 fe_cast->addStreachSchurMatrix(S, aoS);
892 else
895 };
896
900 for (auto &fe : list)
901 if (auto fe_cast = dynamic_cast<EpElementBase *>(fe.second))
902 fe_cast->addBubbleSchurMatrix(S, aoS);
903 else
906 };
907
911 for (auto &fe : list)
912 if (auto fe_cast = dynamic_cast<EpElementBase *>(fe.getSharedPtr().get()))
913 fe_cast->addBubbleSchurMatrix(S, aoS);
914 else
917 };
918
922 for (auto &fe : list)
923 if (auto fe_cast = dynamic_cast<EpElementBase *>(fe.second))
924 fe_cast->addSpatialDispStressSchurMatrix(S, aoS);
925 else
928 };
929
933 for (auto &fe : list)
934 if (auto fe_cast = dynamic_cast<EpElementBase *>(fe.getSharedPtr().get()))
935 fe_cast->addSpatialDispStressSchurMatrix(S, aoS);
936 else
939 };
940
944 for (auto &fe : list)
945 if (auto fe_cast = dynamic_cast<EpElementBase *>(fe.second))
946 fe_cast->addOmegaSchurMatrix(S, aoS);
947 else
950 };
951
955 for (auto &fe : list)
956 if (auto fe_cast = dynamic_cast<EpElementBase *>(fe.getSharedPtr().get()))
957 fe_cast->addOmegaSchurMatrix(S, ao);
958 else
961 };
962
969
973
981 aoSBubble);
983 aoSBubble);
985 aoSBubble);
986
994 aoSOmega);
996 aoSOmega);
998 aoSOmega);
999
1002 &schur_spatial_disp_prb_ptr);
1003 if (auto spatial_disp_data =
1009 aoSw);
1011 aoSw);
1013 aoSw);
1014 } else
1016 "Problem does not have sub-problem data");
1017
1018 } else
1020 "Problem does not have sub-problem data");
1021
1022 } else
1024 "Problem does not have sub-problem data");
1025
1026 } else
1028 "Problem does not have sub-problem data");
1029
1031
1034 using VolOp = VolumeElementForcesAndSourcesCore::UserDataOperator;
1036
1037 EshelbianCore &eP;
1038 boost::shared_ptr<SetPtsData> dataFieldEval;
1039 boost::shared_ptr<VolEle> volPostProcEnergy;
1040 boost::shared_ptr<double> gEnergy;
1041
1043 : eP(ep),
1046 volPostProcEnergy(new
VolEle(ep.mField)), gEnergy(new
double) {
1048 dataFieldEval, "EP");
1049 CHKERRABORT(PETSC_COMM_WORLD,
ierr);
1050
1051 auto no_rule = [](
int,
int,
int) {
return -1; };
1052
1053 auto set_element_for_field_eval = [&]() {
1054 boost::shared_ptr<Ele> vol_ele(dataFieldEval->feMethodPtr.lock());
1055 vol_ele->getRuleHook = no_rule;
1056 vol_ele->getUserPolynomialBase() =
1057 boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase());
1058 vol_ele->getOpPtrVector().push_back(new OpL2Transform());
1059
1060 vol_ele->getOpPtrVector().push_back(
1062 eP.piolaStress, eP.dataAtPts->getApproxPAtPts()));
1063 vol_ele->getOpPtrVector().push_back(
1065 eP.bubbleField, eP.dataAtPts->getApproxPAtPts(), MBMAXTYPE));
1066 vol_ele->getOpPtrVector().push_back(
1068 eP.streachTensor, eP.dataAtPts->getLogStreachTensorAtPts(),
1069 MBTET));
1071 eP.rotAxis, eP.dataAtPts->getRotAxisAtPts(), MBTET));
1072 vol_ele->getOpPtrVector().push_back(
1074 eP.materialGradient, eP.dataAtPts->getBigGAtPts(), MBTET));
1076 eP.spatialDisp, eP.dataAtPts->getSmallWAtPts(), MBTET));
1077 vol_ele->getOpPtrVector().push_back(
1078 new OpCalculateRotationAndSpatialGradient(eP.rotAxis,
1079 eP.dataAtPts));
1080 };
1081
1082 auto set_element_for_post_process = [&]() {
1083 volPostProcEnergy->getRuleHook =
VolRule();
1084 volPostProcEnergy->getUserPolynomialBase() =
1085 boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase());
1086 volPostProcEnergy->getOpPtrVector().push_back(new OpL2Transform());
1087
1088 volPostProcEnergy->getOpPtrVector().push_back(
1090 eP.piolaStress, eP.dataAtPts->getApproxPAtPts()));
1091 volPostProcEnergy->getOpPtrVector().push_back(
1093 eP.bubbleField, eP.dataAtPts->getApproxPAtPts(), MBMAXTYPE));
1094 volPostProcEnergy->getOpPtrVector().push_back(
1096 eP.streachTensor, eP.dataAtPts->getLogStreachTensorAtPts(),
1097 MBTET));
1098 volPostProcEnergy->getOpPtrVector().push_back(
1100 eP.rotAxis, eP.dataAtPts->getRotAxisAtPts(), MBTET));
1101 volPostProcEnergy->getOpPtrVector().push_back(
1103 eP.materialGradient, eP.dataAtPts->getBigGAtPts(), MBTET));
1104 volPostProcEnergy->getOpPtrVector().push_back(
1106 eP.spatialDisp, eP.dataAtPts->getSmallWAtPts(), MBTET));
1107 volPostProcEnergy->getOpPtrVector().push_back(
1108 new OpCalculateRotationAndSpatialGradient(eP.rotAxis,
1109 eP.dataAtPts));
1110 volPostProcEnergy->getOpPtrVector().push_back(
1111 new OpCalculateStrainEnergy(eP.spatialDisp, eP.dataAtPts, gEnergy));
1112 };
1113
1114 set_element_for_field_eval();
1115 set_element_for_post_process();
1116 }
1117
1119
1121
1124
1125
1126
1127
1128
1129
1130
1131
1132
1133 auto get_step = [](auto ts_step) {
1134 std::ostringstream ss;
1135 ss << boost::str(boost::format("%d") % static_cast<int>(ts_step));
1136 std::string s = ss.str();
1137 return s;
1138 };
1139
1140 PetscViewer viewer;
1141 CHKERR PetscViewerBinaryOpen(
1142 PETSC_COMM_WORLD, ("restart_" + get_step(ts_step) + ".dat").c_str(),
1143 FILE_MODE_WRITE, &viewer);
1144 CHKERR VecView(ts_u, viewer);
1145 CHKERR PetscViewerDestroy(&viewer);
1146
1147 CHKERR eP.postProcessResults(1,
"out_sol_elastic_" + get_step(ts_step) +
1148 ".h5m");
1149
1150
1151 *gEnergy = 0;
1152 CHKERR eP.mField.loop_finite_elements(problemPtr->getName(),
"EP",
1153 *volPostProcEnergy);
1154
1155 double body_energy;
1156 MPI_Allreduce(gEnergy.get(), &body_energy, 1, MPI_DOUBLE, MPI_SUM,
1157 eP.mField.get_comm());
1158 MOFEM_LOG_C(
"EP", Sev::inform,
"Step %d time %3.4g strain energy %3.6e",
1159 ts_step, ts_t, body_energy);
1160
1161 auto post_proc_at_points = [&](std::array<double, 3> point,
1162 std::string str) {
1164
1165 dataFieldEval->setEvalPoints(point.data(), point.size() / 3);
1166
1167 struct OpPrint :
public VolOp {
1168
1169 EshelbianCore &eP;
1170 std::array<double, 3> point;
1171 std::string str;
1172
1173 OpPrint(EshelbianCore &ep, std::array<double, 3> &point,
1174 std::string &str)
1175 :
VolOp(ep.spatialDisp,
VolOp::OPROW), eP(ep), point(point),
1176 str(str) {}
1177
1179 DataForcesAndSourcesCore::EntData &data) {
1181 if (type == MBTET) {
1182 if (getGaussPts().size2()) {
1183
1184 auto t_h = getFTensor2FromMat<3, 3>(eP.dataAtPts->hAtPts);
1185 auto t_approx_P =
1186 getFTensor2FromMat<3, 3>(eP.dataAtPts->approxPAtPts);
1187
1193 t_cauchy(
i,
j) = (1. / jac) * (t_approx_P(
i,
k) * t_h(
j,
k));
1194
1195 auto add = [&]() {
1196 std::ostringstream s;
1197 s << str << " elem " << getFEEntityHandle() << " ";
1198 return s.str();
1199 };
1200
1201 auto print_tensor = [](
auto &
t) {
1202 std::ostringstream s;
1204 return s.str();
1205 };
1206
1207 std::ostringstream print;
1209 << add() << "comm rank " << eP.mField.get_comm_rank();
1213 << add() << "coords at gauss pts " << getCoordsAtGaussPts();
1215 << add() << "w " << eP.dataAtPts->wAtPts;
1217 << add() << "Piola " << eP.dataAtPts->approxPAtPts;
1219 << add() << "Cauchy " << print_tensor(t_cauchy);
1220 }
1221 }
1223 }
1224 };
1225
1226 if (auto fe_ptr = dataFieldEval->feMethodPtr.lock()) {
1227
1228 fe_ptr->getOpPtrVector().push_back(new OpPrint(eP, point, str));
1230 ->evalFEAtThePoint3D(
1231 point.data(), 1e-12, problemPtr->getName(), "EP",
1232 dataFieldEval, eP.mField.get_comm_rank(),
1234 fe_ptr->getOpPtrVector().pop_back();
1235 }
1236
1238 };
1239
1240
1241 std::array<double, 3> pointA = {48., 60., 5.};
1242 CHKERR post_proc_at_points(pointA,
"Point A");
1244
1245 std::array<double, 3> pointB = {48. / 2., 44. + (60. - 44.) / 2., 0.};
1246 CHKERR post_proc_at_points(pointB,
"Point B");
1248
1249 std::array<double, 3> pointC = {48. / 2., (44. - 0.) / 2., 0.};
1250 CHKERR post_proc_at_points(pointC,
"Point C");
1252
1254 }
1255 };
1256
1257 boost::shared_ptr<FEMethod> monitor_ptr(
new Monitor(*
this));
1260
1261 CHKERR TSAppendOptionsPrefix(ts,
"elastic_");
1262 CHKERR TSSetFromOptions(ts);
1264}
#define MOFEM_LOG_SYNCHRONISE(comm)
Synchronise "SYNC" channel.
#define MOFEM_LOG_C(channel, severity, format,...)
FTensor::Index< 'm', SPACE_DIM > m
PetscErrorCode DMMoFEMGetProblemPtr(DM dm, const MoFEM::Problem **problem_ptr)
Get pointer to problem data structure.
PetscErrorCode DMCreateMatrix_MoFEM(DM dm, Mat *M)
PetscErrorCode DMMoFEMGetTsCtx(DM dm, MoFEM::TsCtx **ts_ctx)
get MoFEM::TsCtx data structure
#define MOFEM_LOG(channel, severity)
Log.
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
auto getVectorAdaptor(T1 ptr, const size_t n)
Get Vector adaptor.
PetscErrorCode TsMonitorSet(TS ts, PetscInt step, PetscReal t, Vec u, void *ctx)
Set monitor for TS solver.
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
constexpr double t
plate stiffness
structure for User Loop Methods on finite elements
Field evaluator interface.
structure to get information form mofem into EntitiesFieldData
keeps basic data about problem
boost::shared_ptr< SubProblemData > subProblemData
intrusive_ptr for managing petsc objects
FEMethodsSequence loopsIJacobian
FEMethodsSequence & getLoopsMonitor()
Get the loops to do Monitor object.
BasicMethodsSequence preProcessIJacobian
BasicMethodsSequence postProcessIJacobian
Volume finite element base.