8 using Ele = ForcesAndSourcesCore;
9 using VolEle = VolumeElementForcesAndSourcesCore;
12 using PtsHashMap = std::map<std::string, std::array<double, 3>>;
21 :
eP(ep),
dataFieldEval(ep.mField.getInterface<FieldEvaluatorInterface>()
24 ierr = ep.mField.getInterface<FieldEvaluatorInterface>()->buildTree3D(
26 CHKERRABORT(PETSC_COMM_WORLD,
ierr);
28 auto no_rule = [](
int,
int,
int) {
return -1; };
30 auto set_element_for_field_eval = [&]() {
32 boost::shared_ptr<Ele> vol_ele(
dataFieldEval->feMethodPtr.lock());
33 vol_ele->getRuleHook = no_rule;
34 vol_ele->getUserPolynomialBase() =
35 boost::shared_ptr<BaseFunction>(
new CGGUserPolynomialBase());
36 AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(vol_ele->getOpPtrVector(),
39 vol_ele->getOpPtrVector().push_back(
new OpCalculateHVecTensorField<3, 3>(
40 eP.piolaStress,
eP.dataAtPts->getApproxPAtPts()));
41 vol_ele->getOpPtrVector().push_back(
42 new OpCalculateHTensorTensorField<3, 3>(
43 eP.bubbleField,
eP.dataAtPts->getApproxPAtPts(), MBMAXTYPE));
44 vol_ele->getOpPtrVector().push_back(
45 new OpCalculateTensor2SymmetricFieldValues<3>(
46 eP.stretchTensor,
eP.dataAtPts->getLogStretchTensorAtPts(),
48 vol_ele->getOpPtrVector().push_back(
new OpCalculateVectorFieldValues<3>(
49 eP.rotAxis,
eP.dataAtPts->getRotAxisAtPts(), MBTET));
50 vol_ele->getOpPtrVector().push_back(
new OpCalculateVectorFieldValues<3>(
51 eP.spatialL2Disp,
eP.dataAtPts->getSmallWL2AtPts(), MBTET));
54 vol_ele->getOpPtrVector().push_back(
new OpCalculateVectorFieldValues<3>(
55 eP.spatialH1Disp,
eP.dataAtPts->getSmallWH1AtPts()));
56 vol_ele->getOpPtrVector().push_back(
58 eP.spatialH1Disp,
eP.dataAtPts->getSmallWGradH1AtPts()));
60 vol_ele->getOpPtrVector().push_back(
61 new OpCalculateRotationAndSpatialGradient(
eP.dataAtPts));
65 auto set_element_for_post_process = [&]() {
70 boost::shared_ptr<BaseFunction>(
new CGGUserPolynomialBase());
71 AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
75 new OpCalculateHVecTensorField<3, 3>(
76 eP.piolaStress,
eP.dataAtPts->getApproxPAtPts()));
78 new OpCalculateHTensorTensorField<3, 3>(
79 eP.bubbleField,
eP.dataAtPts->getApproxPAtPts(), MBMAXTYPE));
81 new OpCalculateTensor2SymmetricFieldValues<3>(
82 eP.stretchTensor,
eP.dataAtPts->getLogStretchTensorAtPts(),
85 new OpCalculateVectorFieldValues<3>(
86 eP.rotAxis,
eP.dataAtPts->getRotAxisAtPts(), MBTET));
88 new OpCalculateVectorFieldValues<3>(
89 eP.spatialL2Disp,
eP.dataAtPts->getSmallWL2AtPts(), MBTET));
93 new OpCalculateVectorFieldValues<3>(
94 eP.spatialH1Disp,
eP.dataAtPts->getSmallWH1AtPts()));
97 eP.spatialH1Disp,
eP.dataAtPts->getSmallWGradH1AtPts()));
100 new OpCalculateRotationAndSpatialGradient(
eP.dataAtPts));
102 new OpCalculateStrainEnergy(
eP.spatialL2Disp,
eP.dataAtPts,
gEnergy));
106 auto reads_post_proc_data =
112 if (!file.is_open()) {
118 while (std::getline(file, line)) {
119 std::istringstream iss(line);
121 double col2, col3, col4;
123 if (iss >> col1 >> col2 >> col3 >> col4) {
124 MOFEM_LOG(
"EP", Sev::verbose) <<
"Read: " << col1 <<
", " << col2
125 <<
", " << col3 <<
", " << col4;
126 pts_hash_map[col1] = {col2, col3, col4};
128 MOFEM_LOG(
"EP", Sev::error) <<
"Error parsing line: " << line;
139 PetscBool test_cook_flg = PETSC_FALSE;
141 &test_cook_flg, PETSC_NULL),
142 "get post proc points");
145 ptsHashMap[
"Point B"] = {48. / 2., 44. + (60. - 44.) / 2., 0.};
146 ptsHashMap[
"Point C"] = {48. / 2., (44. - 0.) / 2., 0.};
158 MOFEM_LOG(
"EP", Sev::inform) <<
"Monitor postProcess";
160 auto get_step = [](
auto ts_step) {
161 std::ostringstream ss;
162 ss << boost::str(boost::format(
"%d") %
static_cast<int>(ts_step));
163 std::string s = ss.str();
174 CHKERR eP.postProcessResults(1,
"out_sol_elastic_" + get_step(ts_step) +
179 CHKERR eP.mField.loop_finite_elements(problemPtr->getName(),
"EP",
182 double body_energy = 0;
183 MPI_Allreduce(
gEnergy.get(), &body_energy, 1, MPI_DOUBLE, MPI_SUM,
184 eP.mField.get_comm());
185 MOFEM_LOG_C(
"EP", Sev::inform,
"Step %d time %3.4g strain energy %3.6e",
186 ts_step, ts_t, body_energy);
188 auto post_proc_at_points = [&](std::array<double, 3> point,
192 dataFieldEval->setEvalPoints(point.data(), point.size() / 3);
194 struct OpPrint :
public VolOp {
197 std::array<double, 3> point;
209 if (getGaussPts().size2()) {
211 auto t_h = getFTensor2FromMat<3, 3>(
eP.dataAtPts->hAtPts);
213 getFTensor2FromMat<3, 3>(
eP.dataAtPts->approxPAtPts);
220 t_cauchy(
i,
j) = (1. / jac) * (t_approx_P(
i,
k) * t_h(
j,
k));
223 std::ostringstream s;
224 s << str <<
" elem " << getFEEntityHandle() <<
" ";
228 auto print_tensor = [](
auto &
t) {
229 std::ostringstream s;
234 std::ostringstream print;
236 << add() <<
"comm rank " <<
eP.mField.get_comm_rank();
240 << add() <<
"coords at gauss pts " << getCoordsAtGaussPts();
242 << add() <<
"w " <<
eP.dataAtPts->wL2AtPts;
244 << add() <<
"Piola " <<
eP.dataAtPts->approxPAtPts;
246 << add() <<
"Cauchy " << print_tensor(t_cauchy);
255 fe_ptr->getOpPtrVector().push_back(
new OpPrint(
eP, point, str));
256 CHKERR eP.mField.getInterface<FieldEvaluatorInterface>()
257 ->evalFEAtThePoint3D(point.data(), 1e-12, problemPtr->getName(),
261 fe_ptr->getOpPtrVector().pop_back();
269 CHKERR post_proc_at_points(pts.second, pts.first);