159 {
161
162 MOFEM_LOG(
"EP", Sev::inform) <<
"Monitor postProcess";
163
164
165 auto get_ents_on_mesh_skin = [&]() {
166
167 std::map<std::string, Range> boundary_entities_vec;
168
169 auto get_bock_vec_impl = [&](auto block_name) {
171 std::regex(
172
173 (boost::format("%s(.*)") % block_name).str()
174
175 ));
176 };
177
178 auto add_blockset_ents_impl = [&](auto &&vec) {
180 for (auto it : vec) {
181 Range boundary_entities;
183 boundary_entities, true);
184 std::string meshset_name = it->getName();
185 boundary_entities_vec[meshset_name] = boundary_entities;
186 boundary_entities.clear();
187 }
189 };
190
191 for (auto block_name : {"SPATIAL_DISP", "FIX", "CONTACT",
192 "SPATIAL_ROTATION", "NORMAL_DISPLACEMENT"}) {
194 "add blockset entities");
195 }
196
197 return boundary_entities_vec;
198 };
199
200 auto boundary_entities_vec = get_ents_on_mesh_skin();
201
202 std::vector<std::tuple<std::string, Range, std::array<double, 6>>>
203 reactionForces;
204 for (const auto &pair : boundary_entities_vec) {
205 reactionForces.push_back(
206 std::make_tuple(pair.first, pair.second,
207 std::array<double, 6>{0.0, 0.0, 0.0, 0.0, 0.0, 0.0}));
208 }
209
212 };
213 auto face_fe =
214 boost::make_shared<FaceElementForcesAndSourcesCore>(
eP.
mField);
215 auto no_rule = [](
int,
int,
int) {
return -1; };
216 face_fe->getRuleHook = integration_rule_face;
218 EshelbianPlasticity::AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
221
222 auto op_side =
224 face_fe->getOpPtrVector().push_back(op_side);
225 auto side_fe_ptr = op_side->getSideFEPtr();
226 auto base_ptr =
227 boost::make_shared<EshelbianPlasticity::CGGUserPolynomialBase>();
228 side_fe_ptr->getUserPolynomialBase() = base_ptr;
229 CHKERR EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
232 auto piola_scale_ptr = boost::make_shared<double>(1.0);
233 side_fe_ptr->getOpPtrVector().push_back(
236 side_fe_ptr->getOpPtrVector().push_back(
237 new OpCalculateHVecTensorField<SPACE_DIM, SPACE_DIM>(
239 side_fe_ptr->getOpPtrVector().push_back(
240 new OpCalculateHTensorTensorField<3, 3>(
242 side_fe_ptr->getOpPtrVector().push_back(
244 side_fe_ptr->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
246 for (auto &[name, ents, reaction_vec] : reactionForces) {
249 }
250
252 *face_fe);
253
254 for (auto &[name, ents, reaction_vec] : reactionForces) {
255 std::array<double, 6> block_reaction_force{0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
256
257 MPI_Allreduce(reaction_vec.data(), &block_reaction_force, 6, MPI_DOUBLE,
259
260 for (auto &force : block_reaction_force) {
261 if (std::abs(force) < 1e-12) {
262 force = 0.0;
263 }
264 }
266 "EP", Sev::inform,
267 "Step %d time %3.4g Block %s Reaction force [%3.6e, %3.6e, %3.6e]",
268 ts_step, ts_t, name.c_str(), block_reaction_force[0],
269 block_reaction_force[1], block_reaction_force[2]);
271 "Step %d time %3.4g Block %s Moment [%3.6e, %3.6e, %3.6e]",
272 ts_step, ts_t, name.c_str(), block_reaction_force[3],
273 block_reaction_force[4], block_reaction_force[5]);
274 }
275
276 auto get_step = [](auto ts_step) {
277 std::ostringstream ss;
278 ss << boost::str(boost::format("%d") % static_cast<int>(ts_step));
279 std::string s = ss.str();
280 return s;
281 };
282
284 PetscViewer viewer;
285 CHKERR PetscViewerBinaryOpen(
286 PETSC_COMM_WORLD, ("restart_" + get_step(ts_step) + ".dat").c_str(),
287 FILE_MODE_WRITE, &viewer);
288 CHKERR VecView(ts_u, viewer);
289 CHKERR PetscViewerDestroy(&viewer);
290 }
291
293 ".h5m");
294
295
296
297
298 auto get_material_force_tags = [&]() {
300 std::vector<Tag> tag(2);
302 "can't get tag");
304 "can't get tag");
305 return tag;
306 };
308
309
310 bool post_process_skeleton = false;
311#ifndef NDEBUG
312 post_process_skeleton = true;
313#endif
314 if (post_process_skeleton) {
316 1, "out_skeleton_" + get_step(ts_step) + ".h5m", PETSC_NULLPTR,
317 get_material_force_tags());
318 }
319
320
326
327 double body_energy = 0;
328 MPI_Allreduce(
gEnergy.get(), &body_energy, 1, MPI_DOUBLE, MPI_SUM,
331 auto crack_area_ptr = boost::make_shared<double>(0.0);
334 "Step %d time %3.4g strain energy %3.6e crack "
335 "area %3.6e",
336 ts_step, ts_t, body_energy, *crack_area_ptr);
339 if (crack_faces.empty()) {
341 }
344 (boost::format("crack_faces_step_%d.vtk") % ts_step).str(),
345 crack_faces);
348 (boost::format("front_edges_step_%d.vtk") % ts_step).str(),
350 }
351
352 } else {
353 MOFEM_LOG_C(
"EP", Sev::inform,
"Step %d time %3.4g strain energy %3.6e",
354 ts_step, ts_t, body_energy);
355 }
356
357 auto post_proc_at_points = [&](std::array<double, 3> point,
358 std::string str) {
360
361 dataFieldEval->setEvalPoints(point.data(), point.size() / 3);
362
363 struct OpPrint :
public VolOp {
364
366 std::array<double, 3> point;
367 std::string str;
368
370 std::string &str)
371 :
VolOp(ep.spatialL2Disp,
VolOp::OPROW), eP(ep), point(point),
372 str(str) {}
373
375 EntitiesFieldData::EntData &data) {
377 if (type == MBTET) {
378 if (getGaussPts().size2()) {
379
380 auto t_h = getFTensor2FromMat<3, 3>(eP.
dataAtPts->hAtPts);
381 auto t_approx_P =
382 getFTensor2FromMat<3, 3>(eP.
dataAtPts->approxPAtPts);
383
389 t_cauchy(
i,
j) = (1. / jac) * (t_approx_P(
i,
k) * t_h(
j,
k));
390
391 auto add = [&]() {
392 std::ostringstream s;
393 s << str << " elem " << getFEEntityHandle() << " ";
394 return s.str();
395 };
396
397 auto print_tensor = [](
auto &
t) {
398 std::ostringstream s;
400 return s.str();
401 };
402
403 std::ostringstream print;
409 << add() << "coords at gauss pts " << getCoordsAtGaussPts();
411 << add() <<
"w " << eP.
dataAtPts->wL2AtPts;
413 << add() <<
"Piola " << eP.
dataAtPts->approxPAtPts;
415 << add() << "Cauchy " << print_tensor(t_cauchy);
416 }
417 }
419 }
420 };
421
423
424 fe_ptr->getOpPtrVector().push_back(
new OpPrint(
eP, point, str));
426 ->evalFEAtThePoint<SPACE_DIM>(
427 point.data(), 1e-12, problemPtr->getName(),
"EP",
dataFieldEval,
430 fe_ptr->getOpPtrVector().pop_back();
431 }
432
434 };
435
436 auto check_external_strain = [&](std::array<double, 3> point,
439
440 dataFieldEval->setEvalPoints(point.data(), point.size() / 3);
443
446 ->evalFEAtThePoint<SPACE_DIM>(
447 point.data(), 1e-12, problemPtr->getName(),
"EP",
dataFieldEval,
450 }
452 VecSetValue(vec, 0, 0.0, ADD_VALUES);
454 auto add = [&]() {
455 std::ostringstream s;
456 s << str << " elem " << getFEEntityHandle() << " ";
457 return s.str();
458 };
466 VecSetValue(vec, 0, disp_at_point, ADD_VALUES);
467 }
468
470 if (ts_t > 0.0) {
471 CHKERR VecAssemblyBegin(vec);
472 CHKERR VecAssemblyEnd(vec);
473 double error;
474 PetscInt idx = 0;
476 CHKERR VecGetValues(vec, 1, &idx, &error);
477 }
478 MPI_Bcast(&error, 1, MPI_DOUBLE, 0, PETSC_COMM_WORLD);
479
480
481
482 if (std::abs(error - 0.25) > 1e-5) {
484 "Atom test %d failed: wrong displacement.",
atom_test);
485 }
486 }
488 };
489
492 PETSC_NULLPTR);
493
495
498 }
499 } else {
500
502 CHKERR post_proc_at_points(pts.second, pts.first);
504 }
505 }
506
508 }
509
511 }
#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
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
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={})
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