v0.14.0
Loading...
Searching...
No Matches
EshelbianMonitor.cpp
Go to the documentation of this file.
1/** @file
2 @brief Contains definition of EshelbianMonitor class.
3 @ingroup EshelbianPlasticty
4*/
5
6struct EshelbianMonitor : public FEMethod {
7
8 using Ele = ForcesAndSourcesCore;
9 using VolEle = VolumeElementForcesAndSourcesCore;
10 using VolOp = VolumeElementForcesAndSourcesCore::UserDataOperator;
11 using SetPtsData = FieldEvaluatorInterface::SetPtsData;
12 using PtsHashMap = std::map<std::string, std::array<double, 3>>;
13
14 EshelbianCore &eP;
15 boost::shared_ptr<SetPtsData> dataFieldEval;
16 boost::shared_ptr<VolEle> volPostProcEnergy;
17 boost::shared_ptr<double> gEnergy;
19
20 EshelbianMonitor(EshelbianCore &ep)
21 : eP(ep), dataFieldEval(ep.mField.getInterface<FieldEvaluatorInterface>()
22 ->getData<VolEle>()),
23 volPostProcEnergy(new VolEle(ep.mField)), gEnergy(new double) {
24 ierr = ep.mField.getInterface<FieldEvaluatorInterface>()->buildTree3D(
25 dataFieldEval, "EP");
26 CHKERRABORT(PETSC_COMM_WORLD, ierr);
27
28 auto no_rule = [](int, int, int) { return -1; };
29
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(),
37 {HDIV, H1, L2});
38
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(),
47 MBTET));
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));
52
53 // H1 displacements
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()));
59
60 vol_ele->getOpPtrVector().push_back(
61 new OpCalculateRotationAndSpatialGradient(eP.dataAtPts));
63 };
64
65 auto set_element_for_post_process = [&]() {
67 volPostProcEnergy->getRuleHook = VolRule();
68
69 volPostProcEnergy->getUserPolynomialBase() =
70 boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase());
71 AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
72 volPostProcEnergy->getOpPtrVector(), {HDIV, H1, L2});
73
74 volPostProcEnergy->getOpPtrVector().push_back(
75 new OpCalculateHVecTensorField<3, 3>(
76 eP.piolaStress, eP.dataAtPts->getApproxPAtPts()));
77 volPostProcEnergy->getOpPtrVector().push_back(
78 new OpCalculateHTensorTensorField<3, 3>(
79 eP.bubbleField, eP.dataAtPts->getApproxPAtPts(), MBMAXTYPE));
80 volPostProcEnergy->getOpPtrVector().push_back(
81 new OpCalculateTensor2SymmetricFieldValues<3>(
82 eP.stretchTensor, eP.dataAtPts->getLogStretchTensorAtPts(),
83 MBTET));
84 volPostProcEnergy->getOpPtrVector().push_back(
85 new OpCalculateVectorFieldValues<3>(
86 eP.rotAxis, eP.dataAtPts->getRotAxisAtPts(), MBTET));
87 volPostProcEnergy->getOpPtrVector().push_back(
88 new OpCalculateVectorFieldValues<3>(
89 eP.spatialL2Disp, eP.dataAtPts->getSmallWL2AtPts(), MBTET));
90
91 // H1 displacements
92 volPostProcEnergy->getOpPtrVector().push_back(
93 new OpCalculateVectorFieldValues<3>(
94 eP.spatialH1Disp, eP.dataAtPts->getSmallWH1AtPts()));
95 volPostProcEnergy->getOpPtrVector().push_back(
97 eP.spatialH1Disp, eP.dataAtPts->getSmallWGradH1AtPts()));
98
99 volPostProcEnergy->getOpPtrVector().push_back(
100 new OpCalculateRotationAndSpatialGradient(eP.dataAtPts));
101 volPostProcEnergy->getOpPtrVector().push_back(
102 new OpCalculateStrainEnergy(eP.spatialL2Disp, eP.dataAtPts, gEnergy));
104 };
105
106 auto reads_post_proc_data =
107 [](PtsHashMap &pts_hash_map) {
109 std::ifstream file(
110 "points.txt"); // Open the file with the name "data.txt"
111
112 if (!file.is_open()) {
114 }
115
116 std::string line;
117
118 while (std::getline(file, line)) {
119 std::istringstream iss(line);
120 std::string col1;
121 double col2, col3, col4;
122
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};
127 } else {
128 MOFEM_LOG("EP", Sev::error) << "Error parsing line: " << line;
129 }
130 }
131
132 file.close(); // Close the file
134 };
135
136 CHK_THROW_MESSAGE(set_element_for_field_eval(), "set element for field");
137 CHK_THROW_MESSAGE(set_element_for_post_process(), "set element for post");
138
139 PetscBool test_cook_flg = PETSC_FALSE;
140 CHK_THROW_MESSAGE(PetscOptionsGetBool(PETSC_NULL, "", "-test_cook_pts",
141 &test_cook_flg, PETSC_NULL),
142 "get post proc points");
143 if (test_cook_flg) {
144 ptsHashMap["Point A"] = {48., 60., 4.999};
145 ptsHashMap["Point B"] = {48. / 2., 44. + (60. - 44.) / 2., 0.};
146 ptsHashMap["Point C"] = {48. / 2., (44. - 0.) / 2., 0.};
147 }
148 CHK_MOAB_THROW(reads_post_proc_data(ptsHashMap), "read post proc points");
149 }
150
151 MoFEMErrorCode preProcess() { return 0; }
152
153 MoFEMErrorCode operator()() { return 0; }
154
155 MoFEMErrorCode postProcess() {
157
158 auto get_step = [](auto ts_step) {
159 std::ostringstream ss;
160 ss << boost::str(boost::format("%d") % static_cast<int>(ts_step));
161 std::string s = ss.str();
162 return s;
163 };
164
165 PetscViewer viewer;
166 CHKERR PetscViewerBinaryOpen(
167 PETSC_COMM_WORLD, ("restart_" + get_step(ts_step) + ".dat").c_str(),
168 FILE_MODE_WRITE, &viewer);
169 CHKERR VecView(ts_u, viewer);
170 CHKERR PetscViewerDestroy(&viewer);
171
172 CHKERR eP.postProcessResults(1, "out_sol_elastic_" + get_step(ts_step) +
173 ".h5m");
174
175 // Loop boundary elements with traction boundary conditions
176 *gEnergy = 0;
177 CHKERR eP.mField.loop_finite_elements(problemPtr->getName(), "EP",
179
180 double body_energy = 0;
181 MPI_Allreduce(gEnergy.get(), &body_energy, 1, MPI_DOUBLE, MPI_SUM,
182 eP.mField.get_comm());
183 MOFEM_LOG_C("EP", Sev::inform, "Step %d time %3.4g strain energy %3.6e",
184 ts_step, ts_t, body_energy);
185
186 auto post_proc_at_points = [&](std::array<double, 3> point,
187 std::string str) {
189
190 dataFieldEval->setEvalPoints(point.data(), point.size() / 3);
191
192 struct OpPrint : public VolOp {
193
194 EshelbianCore &eP;
195 std::array<double, 3> point;
196 std::string str;
197
198 OpPrint(EshelbianCore &ep, std::array<double, 3> &point,
199 std::string &str)
200 : VolOp(ep.spatialL2Disp, VolOp::OPROW), eP(ep), point(point),
201 str(str) {}
202
203 MoFEMErrorCode doWork(int side, EntityType type,
204 EntitiesFieldData::EntData &data) {
206 if (type == MBTET) {
207 if (getGaussPts().size2()) {
208
209 auto t_h = getFTensor2FromMat<3, 3>(eP.dataAtPts->hAtPts);
210 auto t_approx_P =
211 getFTensor2FromMat<3, 3>(eP.dataAtPts->approxPAtPts);
212
216 const double jac = determinantTensor3by3(t_h);
218 t_cauchy(i, j) = (1. / jac) * (t_approx_P(i, k) * t_h(j, k));
219
220 auto add = [&]() {
221 std::ostringstream s;
222 s << str << " elem " << getFEEntityHandle() << " ";
223 return s.str();
224 };
225
226 auto print_tensor = [](auto &t) {
227 std::ostringstream s;
228 s << t;
229 return s.str();
230 };
231
232 std::ostringstream print;
233 MOFEM_LOG("EPSYNC", Sev::inform)
234 << add() << "comm rank " << eP.mField.get_comm_rank();
235 MOFEM_LOG("EPSYNC", Sev::inform)
236 << add() << "point " << getVectorAdaptor(point.data(), 3);
237 MOFEM_LOG("EPSYNC", Sev::inform)
238 << add() << "coords at gauss pts " << getCoordsAtGaussPts();
239 MOFEM_LOG("EPSYNC", Sev::inform)
240 << add() << "w " << eP.dataAtPts->wL2AtPts;
241 MOFEM_LOG("EPSYNC", Sev::inform)
242 << add() << "Piola " << eP.dataAtPts->approxPAtPts;
243 MOFEM_LOG("EPSYNC", Sev::inform)
244 << add() << "Cauchy " << print_tensor(t_cauchy);
245 }
246 }
248 }
249 };
250
251 if (auto fe_ptr = dataFieldEval->feMethodPtr.lock()) {
252
253 fe_ptr->getOpPtrVector().push_back(new OpPrint(eP, point, str));
254 CHKERR eP.mField.getInterface<FieldEvaluatorInterface>()
255 ->evalFEAtThePoint3D(point.data(), 1e-12, problemPtr->getName(),
256 "EP", dataFieldEval, eP.mField.get_comm_rank(),
257 eP.mField.get_comm_rank(), nullptr, MF_EXIST,
258 QUIET);
259 fe_ptr->getOpPtrVector().pop_back();
260 }
261
263 };
264
265 // Points for Cook beam
266 for (auto &pts : ptsHashMap) {
267 CHKERR post_proc_at_points(pts.second, pts.first);
268 MOFEM_LOG_SEVERITY_SYNC(eP.mField.get_comm(), Sev::inform);
269 }
270
272 }
273};
#define MOFEM_LOG_SEVERITY_SYNC(comm, severity)
Synchronise "SYNC" on curtain severity level.
Definition: LogManager.hpp:352
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:311
static PetscErrorCode ierr
@ QUIET
Definition: definitions.h:208
@ MF_EXIST
Definition: definitions.h:100
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
Definition: definitions.h:595
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
#define CHK_MOAB_THROW(err, msg)
Check error code of MoAB function and throw MoFEM exception.
Definition: definitions.h:576
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
constexpr double t
plate stiffness
Definition: plate.cpp:59
boost::shared_ptr< double > gEnergy
MoFEMErrorCode preProcess()
boost::shared_ptr< VolEle > volPostProcEnergy
MoFEMErrorCode operator()()
boost::shared_ptr< SetPtsData > dataFieldEval
EshelbianMonitor(EshelbianCore &ep)
EshelbianCore & eP
ForcesAndSourcesCore Ele
std::map< std::string, std::array< double, 3 > > PtsHashMap
VolumeElementForcesAndSourcesCore VolEle
VolumeElementForcesAndSourcesCore::UserDataOperator VolOp
MoFEMErrorCode postProcess()
FieldEvaluatorInterface::SetPtsData SetPtsData
@ OPROW
operator doWork function is executed on FE rows
Set integration rule.
VolEle::UserDataOperator VolOp