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 Range boundary_entities;
170 std::string meshset_name = it->getName();
171 std::string block_name_disp = "SPATIAL_DISP";
172 std::string block_name_fix = "FIX";
173 if (meshset_name.compare(0, block_name_disp.size(), block_name_disp) == 0 ||
174 meshset_name.compare(0, block_name_fix.size(), block_name_fix) == 0) {
176 boundary_entities, true);
177
178 boundary_entities_vec[meshset_name] = boundary_entities;
179 boundary_entities.clear();
180 }
181 }
182 return boundary_entities_vec;
183 };
184
185 auto boundary_entities_vec = get_ents_on_mesh_skin();
186
187 std::vector<std::tuple<std::string, Range, std::array<double, 6>>>
188 reactionForces;
189 for (const auto &pair : boundary_entities_vec) {
190 reactionForces.push_back(
191 std::make_tuple(pair.first, pair.second,
192 std::array<double, 6>{0.0, 0.0, 0.0, 0.0, 0.0, 0.0}));
193 }
194
195 auto integration_rule_face = [](int, int,
int approx_order) {
197 };
198 auto face_fe =
199 boost::make_shared<FaceElementForcesAndSourcesCore>(
eP.
mField);
200 auto no_rule = [](int, int, int) { return -1; };
201 face_fe->getRuleHook = integration_rule_face;
206
207 auto op_side =
209 face_fe->getOpPtrVector().push_back(op_side);
210 auto side_fe_ptr = op_side->getSideFEPtr();
211 auto base_ptr =
212 boost::make_shared<EshelbianPlasticity::CGGUserPolynomialBase>();
213 side_fe_ptr->getUserPolynomialBase() = base_ptr;
217 auto piola_scale_ptr = boost::make_shared<double>(1.0);
218 side_fe_ptr->getOpPtrVector().push_back(
221 side_fe_ptr->getOpPtrVector().push_back(
222 new OpCalculateHVecTensorField<SPACE_DIM, SPACE_DIM>(
224 side_fe_ptr->getOpPtrVector().push_back(
225 new OpCalculateHTensorTensorField<3, 3>(
227 side_fe_ptr->getOpPtrVector().push_back(
229 side_fe_ptr->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
231 for (auto &[name, ents, reaction_vec] : reactionForces) {
234 }
235
237 *face_fe);
238
239 for (auto &[name, ents, reaction_vec] : reactionForces) {
240 std::array<double, 6> block_reaction_force{0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
241
242 MPI_Allreduce(reaction_vec.data(), &block_reaction_force, 6, MPI_DOUBLE,
244
245 for (auto &force : block_reaction_force) {
246 if (std::abs(force) < 1e-12) {
247 force = 0.0;
248 }
249 }
251 "EP", Sev::inform,
252 "Step %d time %3.4g Block %s Reaction force [%3.6e, %3.6e, %3.6e]",
253 ts_step, ts_t, name.c_str(), block_reaction_force[0],
254 block_reaction_force[1], block_reaction_force[2]);
256 "Step %d time %3.4g Block %s Moment [%3.6e, %3.6e, %3.6e]",
257 ts_step, ts_t, name.c_str(), block_reaction_force[3],
258 block_reaction_force[4], block_reaction_force[5]);
259 }
260
261 auto get_step = [](auto ts_step) {
262 std::ostringstream ss;
263 ss << boost::str(boost::format("%d") % static_cast<int>(ts_step));
264 std::string s = ss.str();
265 return s;
266 };
267
269 PetscViewer viewer;
270 CHKERR PetscViewerBinaryOpen(
271 PETSC_COMM_WORLD, ("restart_" + get_step(ts_step) + ".dat").c_str(),
272 FILE_MODE_WRITE, &viewer);
273 CHKERR VecView(ts_u, viewer);
274 CHKERR PetscViewerDestroy(&viewer);
275 }
276
278 ".h5m");
279
280
281 bool post_process_skeleton = false;
282 #ifndef NDEBUG
283 post_process_skeleton = true;
284 #endif
285 if (post_process_skeleton) {
286
287 auto get_material_force_tags = [&]() {
289 std::vector<Tag> tag(2);
291 "can't get tag");
293 "can't get tag");
294 return tag;
295 };
296
299 1, "out_skeleton_" + get_step(ts_step) + ".h5m", PETSC_NULLPTR,
300 get_material_force_tags());
301 }
302
303
309
310 double body_energy = 0;
311 MPI_Allreduce(
gEnergy.get(), &body_energy, 1, MPI_DOUBLE, MPI_SUM,
314 auto crack_area_ptr = boost::make_shared<double>(0.0);
317 "Step %d time %3.4g strain energy %3.6e crack "
318 "area %3.6e",
319 ts_step, ts_t, body_energy, *crack_area_ptr);
322 if (crack_faces.empty()) {
324 }
327 (boost::format("crack_faces_step_%d.vtk") % ts_step).str(),
328 crack_faces);
331 (boost::format("front_edges_step_%d.vtk") % ts_step).str(),
333 }
334
335 } else {
336 MOFEM_LOG_C(
"EP", Sev::inform,
"Step %d time %3.4g strain energy %3.6e",
337 ts_step, ts_t, body_energy);
338 }
339
340 auto post_proc_at_points = [&](std::array<double, 3> point,
341 std::string str) {
343
344 dataFieldEval->setEvalPoints(point.data(), point.size() / 3);
345
346 struct OpPrint :
public VolOp {
347
349 std::array<double, 3> point;
350 std::string str;
351
353 std::string &str)
355 str(str) {}
356
358 EntitiesFieldData::EntData &data) {
360 if (type == MBTET) {
362
364 auto t_approx_P =
366
372 t_cauchy(
i,
j) = (1. / jac) * (t_approx_P(
i,
k) * t_h(
j,
k));
373
374 auto add = [&]() {
375 std::ostringstream s;
377 return s.str();
378 };
379
380 auto print_tensor = [](
auto &
t) {
381 std::ostringstream s;
383 return s.str();
384 };
385
386 std::ostringstream print;
394 << add() <<
"w " << eP.
dataAtPts->wL2AtPts;
396 << add() <<
"Piola " << eP.
dataAtPts->approxPAtPts;
398 << add() << "Cauchy " << print_tensor(t_cauchy);
399 }
400 }
402 }
403 };
404
406
407 fe_ptr->getOpPtrVector().push_back(
new OpPrint(
eP, point, str));
409 ->evalFEAtThePoint<SPACE_DIM>(
410 point.data(), 1e-12, problemPtr->getName(),
"EP",
dataFieldEval,
413 fe_ptr->getOpPtrVector().pop_back();
414 }
415
417 };
418
419 auto check_external_strain = [&](std::array<double, 3> point,
422
423 dataFieldEval->setEvalPoints(point.data(), point.size() / 3);
426
429 ->evalFEAtThePoint<SPACE_DIM>(
430 point.data(), 1e-12, problemPtr->getName(),
"EP",
dataFieldEval,
433 }
435 VecSetValue(vec, 0, 0.0, ADD_VALUES);
437 auto add = [&]() {
438 std::ostringstream s;
439 s << str << " elem " << getFEEntityHandle() << " ";
440 return s.str();
441 };
449 VecSetValue(vec, 0, disp_at_point, ADD_VALUES);
450 }
451
453 if (ts_t > 0.0) {
454 CHKERR VecAssemblyBegin(vec);
455 CHKERR VecAssemblyEnd(vec);
456 double error;
457 PetscInt idx = 0;
459 CHKERR VecGetValues(vec, 1, &idx, &error);
460 }
461 MPI_Bcast(&error, 1, MPI_DOUBLE, 0, PETSC_COMM_WORLD);
462
463
464
465 if (std::abs(error - 0.25) > 1e-5) {
467 "Atom test %d failed: wrong displacement.",
atom_test);
468 }
469 }
471 };
472
475 PETSC_NULLPTR);
476
478
481 }
482 } else {
483
485 CHKERR post_proc_at_points(pts.second, pts.first);
487 }
488 }
489
491 }
492
494 }
#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.
#define _IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet having a particular BC meshset in a moFEM field.
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.
FTensor::Tensor2< FTensor::PackPtr< double *, 1 >, Tensor_Dim1, Tensor_Dim2 > getFTensor2FromMat(MatrixDouble &data)
Get tensor rank 2 (matrix) form data matrix.
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
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
virtual MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Operator for bi-linear form, usually to calculate values on left hand side.
EntityHandle getFEEntityHandle() const
Return finite element entity handle.
MatrixDouble & getCoordsAtGaussPts()
Gauss points and weight, matrix (nb. of points x 3)
@ OPROW
operator doWork function is executed on FE rows
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element