v0.14.0
EshelbianMonitor.cpp
Go to the documentation of this file.
1 /** @file
2  @brief Contains definition of EshelbianMonitor class.
3  @ingroup EshelbianPlasticty
4 */
5 
6 struct 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;
18  PtsHashMap ptsHashMap;
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 = [&]() {
31  MoFEMFunctionBegin;
32  boost::shared_ptr<Ele> vol_ele(dataFieldEval->feMethodPtr.lock());
33  vol_ele->getRuleHook = no_rule;
34  vol_ele->getUserPolynomialBase() =
35  boost::make_shared<CGGUserPolynomialBase>();
36  EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
37  vol_ele->getOpPtrVector(), {HDIV, H1, L2}, eP.materialH1Positions,
38  ep.frontAdjEdges);
39 
40  vol_ele->getOpPtrVector().push_back(new OpCalculateHVecTensorField<3, 3>(
41  eP.piolaStress, eP.dataAtPts->getApproxPAtPts()));
42  vol_ele->getOpPtrVector().push_back(
43  new OpCalculateHTensorTensorField<3, 3>(
44  eP.bubbleField, eP.dataAtPts->getApproxPAtPts(), MBMAXTYPE));
45  if (EshelbianCore::noStretch) {
46  vol_ele->getOpPtrVector().push_back(
47  eP.physicalEquations->returnOpCalculateStretchFromStress(
48  eP.dataAtPts, eP.physicalEquations));
49  } else {
50  vol_ele->getOpPtrVector().push_back(
51  new OpCalculateTensor2SymmetricFieldValues<3>(
52  eP.stretchTensor, eP.dataAtPts->getLogStretchTensorAtPts(),
53  MBTET));
54  }
55  vol_ele->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
56  eP.rotAxis, eP.dataAtPts->getRotAxisAtPts(), MBTET));
57  vol_ele->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
58  eP.spatialL2Disp, eP.dataAtPts->getSmallWL2AtPts(), MBTET));
59 
60  // H1 displacements
61  vol_ele->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
62  eP.spatialH1Disp, eP.dataAtPts->getSmallWH1AtPts()));
63  vol_ele->getOpPtrVector().push_back(
64  new OpCalculateVectorFieldGradient<3, 3>(
65  eP.spatialH1Disp, eP.dataAtPts->getSmallWGradH1AtPts()));
66 
67  vol_ele->getOpPtrVector().push_back(
68  new OpCalculateRotationAndSpatialGradient(eP.dataAtPts));
69  MoFEMFunctionReturn(0);
70  };
71 
72  auto set_element_for_post_process_energy = [&]() {
73  MoFEMFunctionBegin;
74  volPostProcEnergy->getRuleHook = VolRule();
75  auto bubble_cache = boost::make_shared<CGGUserPolynomialBase::CachePhi>(
76  0, 0, MatrixDouble());
77  volPostProcEnergy->getUserPolynomialBase() =
78  boost::make_shared<CGGUserPolynomialBase>(bubble_cache);
79  EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
80  volPostProcEnergy->getOpPtrVector(), {HDIV, H1, L2},
81  eP.materialH1Positions, ep.frontAdjEdges);
82  volPostProcEnergy->getOpPtrVector().push_back(
83  new OpCalculateHVecTensorField<3, 3>(
84  eP.piolaStress, eP.dataAtPts->getApproxPAtPts()));
85  volPostProcEnergy->getOpPtrVector().push_back(
86  new OpCalculateHTensorTensorField<3, 3>(
87  eP.bubbleField, eP.dataAtPts->getApproxPAtPts(), MBMAXTYPE));
88  if (eP.noStretch) {
89  volPostProcEnergy->getOpPtrVector().push_back(
90  eP.physicalEquations->returnOpCalculateStretchFromStress(
91  eP.dataAtPts, eP.physicalEquations));
92  } else {
93  volPostProcEnergy->getOpPtrVector().push_back(
94  new OpCalculateTensor2SymmetricFieldValues<3>(
95  eP.stretchTensor, eP.dataAtPts->getLogStretchTensorAtPts(),
96  MBTET));
97  }
98  volPostProcEnergy->getOpPtrVector().push_back(
99  new OpCalculateVectorFieldValues<3>(
100  eP.rotAxis, eP.dataAtPts->getRotAxisAtPts(), MBTET));
101  volPostProcEnergy->getOpPtrVector().push_back(
102  new OpCalculateRotationAndSpatialGradient(eP.dataAtPts));
103 
104  if (auto op = eP.physicalEquations->returnOpCalculateEnergy(eP.dataAtPts,
105  gEnergy)) {
106  volPostProcEnergy->getOpPtrVector().push_back(op);
107  }
108 
109  MoFEMFunctionReturn(0);
110  };
111 
112  auto reads_post_proc_data =
113  [](PtsHashMap &pts_hash_map) {
114  MoFEMFunctionBegin;
115  std::ifstream file(
116  "points.txt"); // Open the file with the name "data.txt"
117 
118  if (!file.is_open()) {
119  MoFEMFunctionReturnHot(0);
120  }
121 
122  std::string line;
123 
124  while (std::getline(file, line)) {
125  std::istringstream iss(line);
126  std::string col1;
127  double col2, col3, col4;
128 
129  if (iss >> col1 >> col2 >> col3 >> col4) {
130  MOFEM_LOG("EP", Sev::verbose) << "Read: " << col1 << ", " << col2
131  << ", " << col3 << ", " << col4;
132  pts_hash_map[col1] = {col2, col3, col4};
133  } else {
134  MOFEM_LOG("EP", Sev::error) << "Error parsing line: " << line;
135  }
136  }
137 
138  file.close(); // Close the file
139  MoFEMFunctionReturn(0);
140  };
141 
142  CHK_THROW_MESSAGE(set_element_for_field_eval(), "set element for field");
143  CHK_THROW_MESSAGE(set_element_for_post_process_energy(),
144  "set element for post energy");
145 
146  PetscBool test_cook_flg = PETSC_FALSE;
147  CHK_THROW_MESSAGE(PetscOptionsGetBool(PETSC_NULL, "", "-test_cook_pts",
148  &test_cook_flg, PETSC_NULL),
149  "get post proc points");
150  if (test_cook_flg) {
151  ptsHashMap["Point A"] = {48., 60., 4.999};
152  ptsHashMap["Point B"] = {48. / 2., 44. + (60. - 44.) / 2., 0.};
153  ptsHashMap["Point C"] = {48. / 2., (44. - 0.) / 2., 0.};
154  }
155  CHK_MOAB_THROW(reads_post_proc_data(ptsHashMap), "read post proc points");
156  }
157 
158  MoFEMErrorCode preProcess() { return 0; }
159 
160  MoFEMErrorCode operator()() { return 0; }
161 
162  MoFEMErrorCode postProcess() {
163  MoFEMFunctionBegin;
164 
165  MOFEM_LOG("EP", Sev::inform) << "Monitor postProcess";
166 
167  auto get_step = [](auto ts_step) {
168  std::ostringstream ss;
169  ss << boost::str(boost::format("%d") % static_cast<int>(ts_step));
170  std::string s = ss.str();
171  return s;
172  };
173 
174  PetscViewer viewer;
175  CHKERR PetscViewerBinaryOpen(
176  PETSC_COMM_WORLD, ("restart_" + get_step(ts_step) + ".dat").c_str(),
177  FILE_MODE_WRITE, &viewer);
178  CHKERR VecView(ts_u, viewer);
179  CHKERR PetscViewerDestroy(&viewer);
180 
181  CHKERR eP.postProcessResults(1, "out_sol_elastic_" + get_step(ts_step) +
182  ".h5m");
183 
184  // Loop boundary elements with traction boundary conditions
185  *gEnergy = 0;
186  TetPolynomialBase::switchCacheBaseOn<HDIV>({volPostProcEnergy.get()});
187  CHKERR eP.mField.loop_finite_elements(problemPtr->getName(), "EP",
188  *volPostProcEnergy);
189  TetPolynomialBase::switchCacheBaseOff<HDIV>({volPostProcEnergy.get()});
190 
191  double body_energy = 0;
192  MPI_Allreduce(gEnergy.get(), &body_energy, 1, MPI_DOUBLE, MPI_SUM,
193  eP.mField.get_comm());
194  MOFEM_LOG_C("EP", Sev::inform, "Step %d time %3.4g strain energy %3.6e",
195  ts_step, ts_t, body_energy);
196 
197  auto post_proc_at_points = [&](std::array<double, 3> point,
198  std::string str) {
199  MoFEMFunctionBegin;
200 
201  dataFieldEval->setEvalPoints(point.data(), point.size() / 3);
202 
203  struct OpPrint : public VolOp {
204 
205  EshelbianCore &eP;
206  std::array<double, 3> point;
207  std::string str;
208 
209  OpPrint(EshelbianCore &ep, std::array<double, 3> &point,
210  std::string &str)
211  : VolOp(ep.spatialL2Disp, VolOp::OPROW), eP(ep), point(point),
212  str(str) {}
213 
214  MoFEMErrorCode doWork(int side, EntityType type,
215  EntitiesFieldData::EntData &data) {
216  MoFEMFunctionBegin;
217  if (type == MBTET) {
218  if (getGaussPts().size2()) {
219 
220  auto t_h = getFTensor2FromMat<3, 3>(eP.dataAtPts->hAtPts);
221  auto t_approx_P =
222  getFTensor2FromMat<3, 3>(eP.dataAtPts->approxPAtPts);
223 
224  FTensor::Index<'i', 3> i;
225  FTensor::Index<'j', 3> j;
226  FTensor::Index<'k', 3> k;
227  const double jac = determinantTensor3by3(t_h);
228  FTensor::Tensor2<double, 3, 3> t_cauchy;
229  t_cauchy(i, j) = (1. / jac) * (t_approx_P(i, k) * t_h(j, k));
230 
231  auto add = [&]() {
232  std::ostringstream s;
233  s << str << " elem " << getFEEntityHandle() << " ";
234  return s.str();
235  };
236 
237  auto print_tensor = [](auto &t) {
238  std::ostringstream s;
239  s << t;
240  return s.str();
241  };
242 
243  std::ostringstream print;
244  MOFEM_LOG("EPSYNC", Sev::inform)
245  << add() << "comm rank " << eP.mField.get_comm_rank();
246  MOFEM_LOG("EPSYNC", Sev::inform)
247  << add() << "point " << getVectorAdaptor(point.data(), 3);
248  MOFEM_LOG("EPSYNC", Sev::inform)
249  << add() << "coords at gauss pts " << getCoordsAtGaussPts();
250  MOFEM_LOG("EPSYNC", Sev::inform)
251  << add() << "w " << eP.dataAtPts->wL2AtPts;
252  MOFEM_LOG("EPSYNC", Sev::inform)
253  << add() << "Piola " << eP.dataAtPts->approxPAtPts;
254  MOFEM_LOG("EPSYNC", Sev::inform)
255  << add() << "Cauchy " << print_tensor(t_cauchy);
256  }
257  }
258  MoFEMFunctionReturn(0);
259  }
260  };
261 
262  if (auto fe_ptr = dataFieldEval->feMethodPtr.lock()) {
263 
264  fe_ptr->getOpPtrVector().push_back(new OpPrint(eP, point, str));
265  CHKERR eP.mField.getInterface<FieldEvaluatorInterface>()
266  ->evalFEAtThePoint3D(point.data(), 1e-12, problemPtr->getName(),
267  "EP", dataFieldEval, eP.mField.get_comm_rank(),
268  eP.mField.get_comm_rank(), nullptr, MF_EXIST,
269  QUIET);
270  fe_ptr->getOpPtrVector().pop_back();
271  }
272 
273  MoFEMFunctionReturn(0);
274  };
275 
276  // Points for Cook beam
277  for (auto &pts : ptsHashMap) {
278  CHKERR post_proc_at_points(pts.second, pts.first);
279  MOFEM_LOG_SEVERITY_SYNC(eP.mField.get_comm(), Sev::inform);
280  }
281 
282  MoFEMFunctionReturn(0);
283  }
284 };
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
CHK_MOAB_THROW
#define CHK_MOAB_THROW(err, msg)
Check error code of MoAB function and throw MoFEM exception.
Definition: definitions.h:589
MOFEM_LOG_SEVERITY_SYNC
#define MOFEM_LOG_SEVERITY_SYNC(comm, severity)
Synchronise "SYNC" on curtain severity level.
Definition: LogManager.hpp:352
EshelbianMonitor::VolOp
VolumeElementForcesAndSourcesCore::UserDataOperator VolOp
Definition: EshelbianMonitor.cpp:10
EshelbianMonitor::operator()
MoFEMErrorCode operator()()
Definition: EshelbianMonitor.cpp:160
EshelbianMonitor::PtsHashMap
std::map< std::string, std::array< double, 3 > > PtsHashMap
Definition: EshelbianMonitor.cpp:12
CHK_THROW_MESSAGE
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
Definition: definitions.h:609
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
EshelbianMonitor::Ele
ForcesAndSourcesCore Ele
Definition: EshelbianMonitor.cpp:8
EshelbianPlasticity::AddHOOps
Definition: EshelbianPlasticity.hpp:1150
VolRule
Set integration rule to volume elements.
Definition: simple_interface.cpp:88
FTensor::Tensor2< double, 3, 3 >
EshelbianMonitor::postProcess
MoFEMErrorCode postProcess()
Definition: EshelbianMonitor.cpp:162
EshelbianMonitor::ptsHashMap
PtsHashMap ptsHashMap
Definition: EshelbianMonitor.cpp:18
FEMethod
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
EshelbianMonitor::volPostProcEnergy
boost::shared_ptr< VolEle > volPostProcEnergy
Definition: EshelbianMonitor.cpp:16
MoFEM::getVectorAdaptor
auto getVectorAdaptor(T1 ptr, const size_t n)
Get Vector adaptor.
Definition: Templates.hpp:31
MOFEM_LOG_C
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:311
double
convert.type
type
Definition: convert.py:64
SetPtsData
FieldEvaluatorInterface::SetPtsData SetPtsData
Definition: dynamic_first_order_con_law.cpp:37
EshelbianMonitor::gEnergy
boost::shared_ptr< double > gEnergy
Definition: EshelbianMonitor.cpp:17
t
constexpr double t
plate stiffness
Definition: plate.cpp:58
EshelbianMonitor::EshelbianMonitor
EshelbianMonitor(EshelbianCore &ep)
Definition: EshelbianMonitor.cpp:20
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
EshelbianMonitor::eP
EshelbianCore & eP
Definition: EshelbianMonitor.cpp:14
EshelbianMonitor::preProcess
MoFEMErrorCode preProcess()
Definition: EshelbianMonitor.cpp:158
EntData
EntitiesFieldData::EntData EntData
Definition: child_and_parent.cpp:37
FTensor::Index< 'i', 3 >
MoFEM::determinantTensor3by3
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition: Templates.hpp:1540
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
UserDataOperator
ForcesAndSourcesCore::UserDataOperator UserDataOperator
Definition: HookeElement.hpp:75
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
EshelbianMonitor::SetPtsData
FieldEvaluatorInterface::SetPtsData SetPtsData
Definition: EshelbianMonitor.cpp:11
EshelbianMonitor::VolEle
VolumeElementForcesAndSourcesCore VolEle
Definition: EshelbianMonitor.cpp:9
MoFEM::Exceptions::ierr
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
EshelbianMonitor
Definition: EshelbianMonitor.cpp:6
EshelbianMonitor::dataFieldEval
boost::shared_ptr< SetPtsData > dataFieldEval
Definition: EshelbianMonitor.cpp:15
QUIET
@ QUIET
Definition: definitions.h:221
k
FTensor::Index< 'k', 3 > k
Definition: matrix_function.cpp:20
MF_EXIST
@ MF_EXIST
Definition: definitions.h:113
convert.int
int
Definition: convert.py:64
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
OpCalculateVectorFieldGradient
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
MoFEM::ForcesAndSourcesCore::UserDataOperator::OPROW
@ OPROW
operator doWork function is executed on FE rows
Definition: ForcesAndSourcesCore.hpp:567
MoFEM::PetscOptionsGetBool
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
Definition: DeprecatedPetsc.hpp:182
Generated by Doxygen 1.8.17 and hosted at University of Glasgow