181 {
183
184 MOFEM_LOG(
"EP", Sev::inform) <<
"Monitor postProcess";
185
186 auto get_step = [](auto ts_step) {
187 std::ostringstream ss;
188 ss << boost::str(boost::format("%d") % static_cast<int>(ts_step));
189 std::string s = ss.str();
190 return s;
191 };
192
193 auto calculate_reaction_forces = [&]() {
195
196
197 auto get_ents_on_mesh_skin = [&]() {
198 std::map<std::string, Range> boundary_entities_vec;
199
200 auto get_block_vec_impl = [&](auto block_name) {
202 std::regex(
203
204 (boost::format("%s(.*)") % block_name).str()
205
206 ));
207 };
208
209 auto add_blockset_ents_impl = [&](auto &&vec) {
211 for (auto it : vec) {
212 Range boundary_entities;
214 boundary_entities, true);
215 std::string meshset_name = it->getName();
216 boundary_entities_vec[meshset_name] = boundary_entities;
217 boundary_entities.clear();
218 }
220 };
221
222 for (auto block_name : {"SPATIAL_DISP", "FIX", "CONTACT",
223 "SPATIAL_ROTATION", "NORMAL_DISPLACEMENT"}) {
225 add_blockset_ents_impl(get_block_vec_impl(block_name)),
226 "add blockset entities");
227 }
228
229 return boundary_entities_vec;
230 };
231
232 auto boundary_entities_vec = get_ents_on_mesh_skin();
233
234 std::vector<std::tuple<std::string, Range, std::array<double, 6>>>
235 reactionForces;
236 for (const auto &pair : boundary_entities_vec) {
237 reactionForces.push_back(std::make_tuple(
238 pair.first, pair.second,
239 std::array<double, 6>{0.0, 0.0, 0.0, 0.0, 0.0, 0.0}));
240 }
241
244 };
245 auto face_fe =
246 boost::make_shared<FaceElementForcesAndSourcesCore>(
eP.
mField);
247 auto no_rule = [](
int,
int,
int) {
return -1; };
248 face_fe->getRuleHook = integration_rule_face;
250 EshelbianPlasticity::AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
253
254 auto op_side =
256 face_fe->getOpPtrVector().push_back(op_side);
257 auto side_fe_ptr = op_side->getSideFEPtr();
258 auto base_ptr =
259 boost::make_shared<EshelbianPlasticity::CGGUserPolynomialBase>();
260 side_fe_ptr->getUserPolynomialBase() = base_ptr;
262 EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
265 auto piola_scale_ptr = boost::make_shared<double>(1.0);
266 side_fe_ptr->getOpPtrVector().push_back(
269 side_fe_ptr->getOpPtrVector().push_back(
270 new OpCalculateHVecTensorField<SPACE_DIM, SPACE_DIM>(
272 piola_scale_ptr));
273 side_fe_ptr->getOpPtrVector().push_back(
274 new OpCalculateHTensorTensorField<3, 3>(
276 side_fe_ptr->getOpPtrVector().push_back(
278 side_fe_ptr->getOpPtrVector().push_back(
279 new OpCalculateVectorFieldValues<3>(
281 for (auto &[name, ents, reaction_vec] : reactionForces) {
284 }
285
288
289 for (auto &[name, ents, reaction_vec] : reactionForces) {
290 std::array<double, 6> block_reaction_force{0.0, 0.0, 0.0,
291 0.0, 0.0, 0.0};
292
293 MPI_Allreduce(reaction_vec.data(), &block_reaction_force, 6, MPI_DOUBLE,
295
296 for (auto &force : block_reaction_force) {
297 if (std::abs(force) < 1e-12) {
298 force = 0.0;
299 }
300 }
302 "Step %d time %3.4g Block %s Reaction force [%3.6e, "
303 "%3.6e, %3.6e]",
304 ts_step, ts_t, name.c_str(), block_reaction_force[0],
305 block_reaction_force[1], block_reaction_force[2]);
307 "Step %d time %3.4g Block %s Moment [%3.6e, %3.6e, %3.6e]",
308 ts_step, ts_t, name.c_str(), block_reaction_force[3],
309 block_reaction_force[4], block_reaction_force[5]);
310 }
312 };
313
314 auto save_restart_file = [&]() {
316
318 PetscViewer viewer;
319 CHKERR PetscViewerBinaryOpen(
320 PETSC_COMM_WORLD, ("restart_" + get_step(ts_step) + ".dat").c_str(),
321 FILE_MODE_WRITE, &viewer);
322 CHKERR PetscViewerBinarySetSkipInfo(viewer, PETSC_TRUE);
323 CHKERR VecView(ts_u, viewer);
324 CHKERR PetscViewerDestroy(&viewer);
325 }
327 };
328
329 auto post_process_skin = [&]() {
332 1, "out_sol_elastic_" + get_step(ts_step) + ".h5m", PETSC_NULLPTR,
333 PETSC_NULLPTR, {}, ts);
335 };
336
337 auto post_process_material_forces = [&]() {
339
342 };
343
344 auto post_process_skeleton_results = [&]() {
346 auto get_material_force_tags = [&]() {
348 std::vector<Tag> tag(2);
350 "can't get tag");
352 "can't get tag");
353 return tag;
354 };
355
356 bool post_process_skeleton = true;
357#ifndef NDEBUG
358 post_process_skeleton = true;
359#endif
360 if (post_process_skeleton) {
362 1, "out_skeleton_" + get_step(ts_step) + ".h5m", PETSC_NULLPTR,
363 get_material_force_tags());
364 }
366 };
367
368 auto calculate_energy = [&]() {
370
376
377 double body_energy = 0;
378 MPI_Allreduce(
gEnergy.get(), &body_energy, 1, MPI_DOUBLE, MPI_SUM,
381 auto crack_area_ptr = boost::make_shared<double>(0.0);
384 "Step %d time %3.4g strain energy %3.6e crack "
385 "area %3.6e",
386 ts_step, ts_t, body_energy, *crack_area_ptr);
389 if (crack_faces.empty()) {
391 }
394 (boost::format("crack_faces_step_%d.vtk") % ts_step).str(),
395 crack_faces);
398 (boost::format("front_edges_step_%d.vtk") % ts_step).str(),
400 }
401
402 } else {
403 MOFEM_LOG_C(
"EP", Sev::inform,
"Step %d time %3.4g strain energy %3.6e",
404 ts_step, ts_t, body_energy);
405 }
407 };
408
409 auto post_proc_at_points = [&](std::array<double, 3> point,
410 std::string str) {
412
413 dataFieldEval->setEvalPoints(point.data(), point.size() / 3);
414
415 struct OpPrint :
public VolOp {
416
418 std::array<double, 3> point;
419 std::string str;
420
422 std::string &str)
423 :
VolOp(ep.spatialL2Disp,
VolOp::OPROW), eP(ep), point(point),
424 str(str) {}
425
427 EntitiesFieldData::EntData &data) {
429 if (type == MBTET) {
430 if (getGaussPts().size2()) {
431
432 auto t_h = getFTensor2FromMat<3, 3>(eP.
dataAtPts->hAtPts);
433 auto t_approx_P =
434 getFTensor2FromMat<3, 3>(eP.
dataAtPts->approxPAtPts);
435
441 t_cauchy(
i,
j) = (1. / jac) * (t_approx_P(
i,
k) * t_h(
j,
k));
442
443 auto add = [&]() {
444 std::ostringstream s;
445 s << str << " elem " << getFEEntityHandle() << " ";
446 return s.str();
447 };
448
449 auto print_tensor = [](
auto &
t) {
450 std::ostringstream s;
452 return s.str();
453 };
454
455 std::ostringstream print;
461 << add() << "coords at gauss pts " << getCoordsAtGaussPts();
463 << add() <<
"w " << eP.
dataAtPts->wL2AtPts;
465 << add() <<
"Piola " << eP.
dataAtPts->approxPAtPts;
467 << add() << "Cauchy " << print_tensor(t_cauchy);
468 }
469 }
471 }
472 };
473
475
476 fe_ptr->getOpPtrVector().push_back(
new OpPrint(
eP, point, str));
478 ->evalFEAtThePoint<SPACE_DIM>(
479 point.data(), 1e-12, problemPtr->getName(),
"EP",
dataFieldEval,
482 fe_ptr->getOpPtrVector().pop_back();
483 }
484
486 };
487
488 auto check_external_strain = [&](std::array<double, 3> point,
491
492 dataFieldEval->setEvalPoints(point.data(), point.size() / 3);
495
498 ->evalFEAtThePoint<SPACE_DIM>(
499 point.data(), 1e-12, problemPtr->getName(),
"EP",
dataFieldEval,
502 }
504 VecSetValue(vec, 0, 0.0, ADD_VALUES);
506 auto add = [&]() {
507 std::ostringstream s;
508 s << str << " elem " << getFEEntityHandle() << " ";
509 return s.str();
510 };
518 VecSetValue(vec, 0, disp_at_point, ADD_VALUES);
519 }
520
522 if (ts_t > 0.0) {
523 CHKERR VecAssemblyBegin(vec);
524 CHKERR VecAssemblyEnd(vec);
525 double error;
526 PetscInt idx = 0;
528 CHKERR VecGetValues(vec, 1, &idx, &error);
529 }
530 MPI_Bcast(&error, 1, MPI_DOUBLE, 0, PETSC_COMM_WORLD);
531
532
533
534 if (std::abs(error - 0.25) > 1e-5) {
536 "Atom test %d failed: wrong displacement.",
atom_test);
537 }
538 }
540 };
541
542 CHKERR save_restart_file();
543 CHKERR calculate_reaction_forces();
544 CHKERR calculate_energy();
545
546 CHKERR post_process_material_forces();
547 CHKERR post_process_skin();
548 CHKERR post_process_skeleton_results();
549
552 PETSC_NULLPTR);
553
555
558 }
559 } else {
560
562 CHKERR post_proc_at_points(pts.second, pts.first);
564 }
565 }
566
568 }
#define MOFEM_LOG_SEVERITY_SYNC(comm, severity)
Synchronise "SYNC" on curtain severity level.
#define MOFEM_LOG_C(channel, severity, format,...)
@ MOFEM_ATOM_TEST_INVALID
virtual MoFEMErrorCode loop_finite_elements(const std::string problem_name, const std::string &fe_name, FEMethod &method, boost::shared_ptr< NumeredEntFiniteElement_multiIndex > fe_ptr=nullptr, MoFEMTypes bh=MF_EXIST, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr(), int verb=DEFAULT_VERBOSITY)=0
Make a loop over finite elements.
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.
auto createVectorMPI(MPI_Comm comm, PetscInt n, PetscInt N)
Create MPI Vector.
MoFEMErrorCode determinantTensor3by3(T1 &t, T2 &det)
Calculate determinant 3 by 3.
constexpr double t
plate stiffness
static constexpr int approx_order
MoFEMErrorCode postProcessResults(const int tag, const std::string file, Vec f_residual=PETSC_NULLPTR, Vec var_vec=PETSC_NULLPTR, std::vector< Tag > tags_to_transfer={}, TS ts=PETSC_NULLPTR)
static PetscBool crackingOn
MoFEMErrorCode calculateFaceMaterialForce(const int tag, TS ts)
MoFEMErrorCode calculateCrackArea(boost::shared_ptr< double > area_ptr)
MoFEMErrorCode postProcessSkeletonResults(const int tag, const std::string file, Vec f_residual=PETSC_NULLPTR, std::vector< Tag > tags_to_transfer={})
const std::string skinElement
boost::shared_ptr< Range > crackFaces
boost::shared_ptr< Range > frontEdges
virtual moab::Interface & get_moab()=0
virtual MPI_Comm & get_comm() const =0
virtual int get_comm_rank() const =0
